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!>■ ; ABSTRACT 

o 

^^ . We study the ionization and thermal evolution of the intergalactic medium during 

[^ , the epoch of He ii reionization by means of radiation hydrodynamical cosmological sim- 

^' I ulations. We post-process baryonic density fields from a standard optically-thin IGM 

simulation with a homogeneous galaxy-dominated UV background (UVB) which reion- 
izes Hi and Hei at z=6.5 but does not have any contribution to the ionization of Heii . 
Therefore, we suppress the Heii photoheating contribution to the gas temperature due 
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1-^ ' to the homogeneous UVB. Quasars are introduced as point sources throughout the 100 

Mpc simulation volume located at cold dark matter (CDM) density peaks consistent 
with the Pei luminosity function. We assume an intrinsic quasar spectrum J(i^) oc v~^'^ 
c/3 , and a luminosity proportional to the halo mass. We evolve the spatial distribution of 

the Hell ionizing radiation field at hv = 4, 8, and 16 Ryd using a time-implicit variable 
tensor Eddington factor radiative transfer scheme. Simultaneously, we also solve for the 
^ ' local ionization of Heii to Heii and the associated photoheating of the gas including 

opacity effects. We find that the percolation of the He in regions is essentially complete 

Q^ . by z=2.5. When comparing to a self-consistent optically thin simulation at the same 

redshift, in which Hell is also ionized by the uniform UVB, we find that inclusion of 
opacity effects results in higher IGM temperature by a factor of approximately 1.7 at 

t~^ I the mean gas density level. We construct synthetic absorption line spectra from which 

we derive statistical parameters of the Heii Lya forest . We use 300 long (Az = 0.2) 
random lines of sight to compute at z = 2.5ib0.1 a mean Heii Lya line transmission 

^ ■ of -F = 0.304 lb 0.002. The error corresponds to a significant one standard deviation in 

j^ ■ the transmitted fiux due to the sightline to sightline variance equal to ~ 11% the mean 

value. The opacity effect on the gas temperature is shown by comparing the broadening 
width of the H i and He ii Lya lines to the results from the self-consistent optically thin 
simulation. We find a shift by approximately 1.25 km/s to higher b-parameter values 
for both H I and He ii . Finally, we estimate the relative broadening width between the 
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two forests and find that the Hell median b-parameter is about 0.8 times the median 
Hi broadening width. This imphes that the Heii absorbers are physically extended 
consistent with conclusions from observed lines of sight. 

Subject headings: cosmology: hydro dynamical simulations: radiation transfer 



1. Introduction 

The thermal evolution of the IGM after reionization is governed chiefly by photoheating (Efs- 
tathiou 1992; Miralda-Escude & Rees 1994). Models for the thermal evolution of the IGM (Miralda- 
Escude &: Ostriker 1990; Giroux & Shapiro 1996; Abel &: Haehnelt 1999) generally assume that 
Hi and Hei are nearly fully ionized by z ~ 6 by star forming galaxies (Fan, Carilli & Keating 
2006), while Heii ionizes somewhat later at z ~ 3 by quasars due to Hen 's higher ionization po- 
tential and recombination rate (Sokasian, Abel & Hernquist 2002; hereafter SAH). Observational 
support for late Heii reionization is summarized in SAH. Abel & Haehnelt (1999) emphasized the 
importance of opacity effects during reionization in establishing the post-reionization temperature 
of the gas. They showed that models using the optically thin expression for the photoheating rate 
during Hell reionization underestimate the IGM temperature at mean density by a factor of ~ 2. 

Hydrodynamical cosmological simulations of the Lya forest (Gen et al. 1994; Zhang, Anninos 
& Norman 1995; Hernquist et al. 1996; Miralda-Escude et al. 1996; Zhang et al. 1997, 1998; 
Theuns et al. 1998) generally adopt the optically thin expression for photoheating for simplicity 
and computational economy. This is a reasonable assumption for the Hi and Hei photoheating 
at z ~ 3 due to the low opacity of the IGM, but not for the He ii photoheating as He ii is in the 
process of being ionized by the percolation of Hem spheres (SAH). These simulations therefore 
underestimate the temperature of the IGM during the epoch of helium reionization. The stan- 
dard approach taken in these simulations is to assume a homogeneous photoionizing background 
which evolves with redshift consistent with observed quasar and galaxy counts, such as that by 
Haardt & Madau (2001). The ionization and thermal state of the baryonic gas is then computed 
self-consistently with its dynamics by solving the equations of hydrodynamic cosmology including 
radiative heating and cooling (Cen 1992; Katz, Weinberg & Hernquist 1996; Anninos et al. 1997). 
It is found that the temperature of the low density IGM is determined almost entirely by photo- 
heating balancing adiabatic cooling due to cosmic expansion. This results in a tight relationship 
between gas temperature and density, the so-called equation of state of the IGM (Hui & Gnedin 
1997): 

T = To/\l^ (1) 

where A = p/p is the gas overdensity and (3 is redshift dependent but in the range < I3{z) < 0.6. 
Abel & Haehnelt (1999) show that opacity effects during Heii reionization raises Tq and reduces 
(3 relative to an optically thin calculation. A technique often employed to include opacity effects 
within hydrodynamic simulations is to simply multiply the Hell photoheating rate by a constant 



factor ^fjg jj > 1, with a value of 2-4 being sufficient to match observed temperatures (e.g., Bryan 
& Machacek 2000; Jena et al. 2005). 

The combination of high spectral resolution quasar absorption line observations and hydro- 
dynamical cosmological simulations provide a means for measuring the thermal evolution of the 
IGM. The thermal state of the gas is generally deduced from H i Lya linewidths (b- values) in high 
resolution spectra (Ranch et al. 1997; Schaye et al. 1999, 2000; Bryan & Machacek 2000; Theuns 
et al. 2000; Bolton et al. 2005; Jena et al. 2005) although the flux power spectrum has also been 
employed (Zaldarriaga, Hui & Tegmark 2001). If the temperature of the IGM alone determined 
the b-values, then it would be straightforward to measure the temperature from high resolution 
spectra. However, Hubble broadening is always of the same order as the thermal broadening in 
Lya forest absorption lines (Schaye 1999), hence the need for comparison with simulations. The- 
uns, Schaye & Haehnelt (2000) used the b-parameter distribution to measure the temperature of 
the IGM at z=3.25, finding Tq > 15, 000 K. Schaye et al. (2000) used the lower cutoff in the 
linewidth-column density scatter diagram to measure the temperature evolution of the IGM over 
the redshift range 2-4.5. They found evidence of late reheating at z~3, which they ascribed to late 
Hen reionization by quasars. Bryan &: Machacek (2000) independently explored the same diagnos- 
tics and found that temperature estimates were sensitive to the assumed cosmology, in particular 
the amplitude of mass fiuctuations on a few Mpc scales. Zaldarriaga, Hui & Tegmark (2001) used 
the falloff of the Lya forest fiux power spectrum at small scales to measure the IGM temperature, 
and found Tq '^ 2 x 10^ K, dependent on their assumed (5. 

A proper calculation of Hell reionization would treat quasars as point sources within the 
simulation volume and solve the equation of radiative transfer for the spatial distribution of the 
ionizing background coupled self-consistently to the dynamical, thermal, and ionization evolution 
of the gas. This was done in an approximate way by SAH, who post-processed a series of density 
fields taken from a SPH hydrodynamical simulation of the IGM using the GADGET code. The 
hydrodynamical simulation used the optically thin prescription for ionizing and heating the gas, 
assuming all ionization states of hydrogen and helium were in ionization equilibrium with the UVB 
of Haardt &; Madau (1996). In the post-processing step, the helium reionization calculation was 
recomputed for an evolving quasar source population treated as point sources within the volume. 
Quasars were assumed to have a constant lifetime of 10 yr. Every 10 yr, peaks within the 
density field of a suitably chosen data dump from the hydrodynamic simulation were populated 
with quasars consistent with an empirical luminosity function. Around each point source the static 
equation of radiative transfer for He ii ionizing photons was solved using a photon-conserving ray- 
casting scheme. The ionization state of Heii was then updated ignoring thermal feedback to the 
gas. SAH found that for reasonable parameter choices, Heii reionization occurred in the range 
3 < z < 4 consistent with observations. 

In this paper we present a simulation which is similar in spirit to SAH, but with several 
important differences. We also postprocess a series of snapshots from a hydrodynamic cosmological 
simulation with a radiative transfer code, however we keep track of the photoheating of the IGM 



as expanding Hem spheres percolate and eventually merge. Our hydrodynamic simulation was 
performed on an Eulerian grid of 512^ cells versus SAH's 224^ in the same volume, giving us ~ 12 
times the mass resolution and ~ 2.3 times the spatial resolution in the low density IGM. We have 
also separated the effects of stellar and quasar populations of the ionization evolution of the IGM 
differently from SAH. Our hydrodynamic simulation computes the ionization of H i and He i due to 
stellar sources only using the homogeneous UVB of Haardt & Madau (2001) arising from galaxies 
"GAL". In the post-processing step, quasars are treated as point sources, which ionize Hell to 
He III . Although our radiative transfer scheme is based on completely different spatial and angular 
discretizations as SAH, the important difference for the purpose of this paper is that we solve the 
RT equation at three frequencies /iz/=4, 8 and 16 Ryd in order to evaluate, albeit crudely, the local 
Hell photoheating rate taking the processed QSO spectrum into account. 

The outline of this paper is as follows. In §2 we describe the cosmological realization for the 
hydrodynamic simulation of homogeneous Hi and He I reionization described in §3. In §2 we also 
introduce some physical concepts relating to late He II reionization and also describe our treatment 
of quasars in the simulation volume. In §4 we describe our method for simulating inhomogeneous 
Hen reionization. In §4.1 our radiative transfer scheme is detailed, in §4.2 we present our single 
species ionization model for Hen , and in §4.3 we present results in terms of globally integrated 
quantities. Since we have split the calculation into two phases, the homogeneous reionization of 
Hi and Hei , followed by the inhomogeneous reionization of Heii , our calculation is not fully 
self-consistent. Although we keep track of the late reheating, we do not modify the underlying 
density fields nor do we alter temperature-dependent recombination rates, which will affect the 
detailed ionization state of the gas. In §5 we present an analysis of our main result, which is the 
late reheating of the IGM due to inhomogeneous Heii reionization, as well as the importance of 
neglecting these coupling effects. We show that these effects, while present, are small, and do not 
seriously undermine our estimate of the reheating. In §6 we present observational signatures of 
Hell reionization based on synthetic Hi and Hell Lya absorption line spectra derived from our 
simulation. In §7 we summarize our main results and conclude. In a series of appendices we derive 
the rate equation for species fraction in an expanding universe (Appendix A), provide more detail 
on the affect of neglecting Heii photoheating on the ionization state of the gas (Appendix B), and 
document the reaction rate coefficients we use (Appendix C). 



2. Simulations 

2.1. Cosmic Realization 

In this work, we present the results from post-processing the redshift evolution of a cosmic 
realization computed with the Eulerian cosmological hydrodynamic code Enzo (Bryan & Norman 
1997; Norman & Bryan 1999; O'Shea et al. 2004). The box of size 67/i~^ Mpc comoving was 
evolved in a flat (il = 1) ACDM cosmology on a unigrid mesh of 512^ grid cells and 512^ dark 



matter particles from z=99 to z ~ 2. We used the initial power spectrum of matter fluctuations 
by Eisenstein &: Hu (1999) to initialize the calculation which was then computed forward under 
the Zel'dovich approximation (Bertschinger & Gelb 1991). Our choice of cosmological parameters 
is (Tg = 0.8, Oa = 0.7, r^ft = 0.04, n^ = 1 with a present day Hubble constant of /i = 0.67 in 
units of 100 km/s/Mpc. With these parameters our cosmic realization has a mesh resolution of 
130h~^ ~ 200 kpc within the 100 Mpc cube, and ~ 1.6h^^ x 10® Mq dark matter particle mass. In 
Figure ([T]) we show a volumetric rendering of the cosmic gas distribution as mapped by the baryon 
overdensity in our simulation at z=2.6. 

One of the fundamental simplifications used here, which places limits on the validity of our 
results is that the ionization of H i and He i are treated entirely differently than He ii . The moti- 
vation is that of numerical simplicity and calculation speed. From the IGM cosmological evolution 
standpoint, hydrogen and neutral helium are believed to be globally ionized at an earlier cosmic 
epoch than He ii . In the simulation discussed here, the photo-ionization/heating rates of H i and 
He I are computed self-consistently during the hydrodynamical calculation in the optically thin 
regime. The premise of our argument is that, if one adopts a picture where neutral hydrogen and 
neutral helium reionization is completed by 2: ~ 6, such species will have large ionization fractions 
by the time He 11 reionization occurs, which theoretical models and observations place at 2: ~ 3 — 2. 

The reionization of the H i /He i species is achieved by the evolving uniform metagalactic flux 
due to stellar sources as computed in Haardt & Madau (2001). The uniform ultra-violet background 
photo-ionizes and photo-heats the IGM, however it is prevented by hand to radiatively alter the 
Hell abundance. We do so by suppressing the ionization/heating rates in Enzo that control the 
Hell <-> Hem chemistry. The latter is computed separately by our inhomogeneous point-source 
distribution of QSOs which is discussed in ^2.31 



The initial calculation represents an undisturbed cosmological ensemble of gas fields ionized 
by a distributed galaxy population. The local QSO component then acts as a perturbation in the 
amplitude of the radiative energy that further ionizes and photoheats the diffuse IGM. The end 
result is that the He 11 reionization proceeds, under the limitations discussed above, via the mergers 
of individual Hem I-fronts during the cosmic epoch that spans the redshift interval z = 6 — 2. 



2.2. Why Late He 11 reionization 

In this section we introduce some definitions and basic results related to He 11 reionization that 
we will refer to later on. The softness parameter of the cosmic radiation is defined as S = t^ "' . In 
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the optical thin limit under a single power law profile for the radiation spectrum, J^ = Jgi2{ , , 
where J912 denotes the volume averaged mean intensity at the hydrogen ionization threshold, we 
can estimate the photoionization rates to be Vui = Fi = "^^h^"^ i+q, a-iid ThbII = r2 = 
^^'■^^Heii ^1 4-Qg_ Dividing the last two relations yields an estimate for the softness parameter 
in the optical thin limit equal to 5" = 4^^°'i . Therefore, the ability of a background radiation 



field to ionize the Hi and Heii species depends on the spectral slope. A soft radiation field would 
have large spectral slope and would primarily favor hydrogen ionization. Hell ionization therefore 
requires a hard ionizing spectrum. 

The available sources of radiation at redshifts z < 10 are of two types. Stellar sources, as- 
sociated with radiation from galaxies forming in the cosmic medium at such redshifts, have large 
softness parameter values of S stellar ~ 4 x 10^. QSOs on the other hand, have much smaller softness 
parameter values of Sqso ~ 50 which makes them ideal for He ii ionization. However, the number 
density of QSOs, as measured by observations sharply rises in the interval 6 > z > 3 and reaches a 
peak at about z ss 3 (Pei 1995). The highest redshift QSO observed to date is at Zem = 6.56 (Hu, 
Cowie &: McMahon 2002). Direct observations of QSO in Lya spectra are very difficult because 
of the hydrogen Gunn- Peterson effect, the optical depth manifestation of hydrogen neutrality at 
high redshifts (Fan, Carilli & Keating 2006). This suggests that either the quasars have a well 
established population by z ss 6 but nevertheless obscured by the large optical depth or that their 
formation began at that epoch. However, when their emissivity is computed at redshifts just past 
the Gunn-Peterson trough, their numbers and spectrum shape is insufficient to produce the hy- 
drogen reionization thought to be primarily achieved by a stellar component in the UVB which 
ionizes hydrogen and neutral helium to the singly ionized state. Because the stellar component 
of the radiation background is not an efficient He ii ionizer. He ii reionization is expected to occur 
later; at times when a sufficient number hard photons becomes available. 

The Gunn-Peterson optical depth of Heii for a uniform IGM can be expressed in the same 
way as the corresponding Hi and that leads to the ratio ''^^"/■^'^ = HM^iLZHeu. = ISjml^ xhe 
ratio of optical depths is known as the R-factor and is used in observations to measure the relative 
properties between the Hell and Hi Lya forest spectra. Observations find that in transmission 
spectra longward of quasars seen at emission redshifts z ~ 3 typically measure values of the R- 
factor between i? ~ 10 — 100 within 6z ~ 0.5 (Reimers et al. 1997; Heap et al. 2000; Kriss et 
al. 2001). This suggests the presence of a HeH Gunn-Peterson trough at redshifts z > 2.5. In the 
optical thin limit the R-factor and the softness parameter relation can be derived. 

In ionization equilibrium the ratio of nHeii/nni can be computed if we assume that the higher 
ionization states of both species (Hem , Hii ) dominate. In that case, for hydrogen we get XiTi = 
neai(l — xi) where xi = nni/nn and ai is Hii the recombination coefficient. Similarly for Heii we 
can also write ^2r2 = '1-602(1 — V'2)) where ^2 = nHell/nue and 02 is the Hell recombination 
coefficient. The above two relations yield Xi/^2 — r^^> for xi ^ 1 and ^2 *^ 1- The last 
equation allows for the determination of the ratio nHi/nueii = ^^n" — 12 x p^"^ . Therefore, for 
— ~ 5 (at T ~ 10^ K) the ratio of number densities becomes SJi^^L ~ jt;S. From the last relation it 
follows that i? ~ 0.1 S. For 5 ~ 4 x 10^ (stellar radiation) the relation predicts R ~ 400. Similarly 
for quasar radiation, with S ~ 50, i? ~ 5. In addition, we can compute that since the quantity 
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neCi{T)/T = 7-«''"y'7-'"ec^ ^]^g ratio of the ionization to recombination time scales, ?i'' = 55" -^is 
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Typical softness parameter values S > 100 then yield "eP > 5 x lO^^r^c. The last relation 



shows that, due to the larger recombination coefficient and the fewer number of photons available for 
ionization, it is more difficult to ionize Hen compared to Hi . The shorter (longer) recombination 
(ionization) time scale effectively restricts Heii reionization to take place at a slower rate. 

Large values in the R-factor is suggestive of large optical depths in He ii and/or small optical 
depths in neutral hydrogen. Comparable optical depths at late redshifts, which correlate to small 
cosmic neutral hydrogen fraction and therefore small cosmic fraction in Hell , yield small values 
for the R-factor. The observations can therefore infer a range in the observed softness parameter 
from measuring the ratio of optical depths in the Hell and Hi line forest. The observed range of 
R=10-100 consequently suggests that 8=100-1000 between z=2-3 although the upper limit is an 
overestimate because we assumed ip2 ■^ ^- Because such observed values are sampled in trans- 
mission spectra that probe the ionization phase transition of Heii to Hem , we can infer that at 
that epoch the galaxy dominated ultraviolet background is gradually being replaced by a quasar 
dominated type. 

In ^ we show that the R-factor evolves with redshift from large to small values as it is com- 
puted in transmission through the computational volume. The quantity that is actually computed 
is the r] parameter defined as the column density ratio between Heii and Hi , r/ = jy"" - The 
R-factor evolution and value follows directly as R = j. The gradual evolution of the R-factor and 
r] parameter from large values is indicative of the cosmic evolution of the softness parameter from 
a stellar to a quasar dominated type. 



2.3. Quasar Placement and Evolution Pre-processing 

For the quasar placement and evolution in the simulation we use the quasar luminosity function 
by Pei et al. (1995), shown in Equation ([2]). This luminosity function was also used by Haardt 
&: Madau (1996,2001) to derive the quasar emissivity and the volume average photoionization and 
heating rates used in the simulation for every species other than Heii . 

'^^'^'^^-[L/L,(z)]ft + [L/L,(z)]ft ^'> 



L,(z) = L,(0)(l + z)"« 



g^Z _^ g^Z, 



Our first requirement is the placement of a single QSO in the computational volume at z=6.5. 
The redshift is for all practical purposes a matter of choice, since there is no accurate prediction of 
when the first QSO appears in the universe. We adopt a scenario where quasars become visible in 
observations after the epoch of hydrogen reionization is completed by z ~ 6.5 due a soft component 
in the ultra-violet background most likely associated with dwarf galaxy formation. In addition, we 
are constrained by the mass resolution of our simulation, which cannot resolve halos smaller than 



~ 5 X 10^ Mq. Therefore, quasar point sources have to be placed in the centers of halos above this 
cutoff. We do so by identifying the evolution of a list of halo centers through the computational 
volume from z = 6.5 to z = 2. The number of quasars per redshift interval is determined by the 
luminosity function in Equation ([2]) and the functional form of the mass-halo to quasar luminosity 
relation. 

We assume an intrinsic quasar spectrum J^ = J9i2i-[r~)~"'' for the LyC part of the spec- 
trum with a spectral index aq = 1.8. Therefore, the number flux of emitted LyC photons, in 
photons / s / cm? , is then hph = ^ ^^ Ji° "^ ^ "^ de. In this relation, e = ^^'^ . The integration 
yields hph = ^^^. Similarly the LyC energy flux is hyc = ^™t^, which yields hyc = 
{hiygi2)hpfi-^ — =^ Li = {hi^9i2) Nph-^ — . In the last equation, Li and Nph represent the total 
LyC luminosity (ergs/s) and emitted photon rate (photons/s) respectively per point source (QSO) 
above the Hi ionization threshold of 1 Ryd (13.6 eV). This allows for a parametrization of the 
emitted ionizing flux based on the number of LyC photons rather than energy. The motivation 
is entirely for consistency with the photon-conserving schemes in simulating reionization (Abel & 
Haehnelt 1999; Sokasian, Abel & Haehnelt 2001, 2002; Ciardi et al. 2003; Whalen, Abel & Norman 
2004). The total luminosity for the Heii ionizing radiation above 4 Ryd (54.4 eV) is then obtained 
through L4 = 4~'"«"'"^) x Li. The luminosity of each source is determined by the dark matter mass 
of the halo that initially creates it. For a dark matter halo of mass M^aio we put in a UV source 
emitting Np^ = 10^^ x ^^'° Hi ionizing photons/s. This is equivalent to placing a Np^ = 10^^ 
ph/ s mini-quasar inside a 10® Mq dark matter halo which is the prescription used in Abel &: 
Haenhelt (1999). Sokasian et al. (2001) investigated an array of QSO placement methods in the 
computational volume and found that the results are largely insensitive to the choice. We adopt 
the linear relation for convenience. In this work, the most massive halo computed at at z=2.5 has 
mass Mfiaio ~ 6 x 10^^ Mq. If a QSO source is placed there then it will emit LyC photons at a rate 
of Nph = 6 X 10^^ s^^. For an input spectrum with slope Oq = 1.8 the photon rate corresponds to a 
LyC luminosity of Li = 2.9 x 10^^ ergs/s and L4 = 6.1 x 10^^ ergs/s for Hi and He 11 respectively. 

The placement of the point sources in the volume is dynamical in nature. The list of dark 
matter halos is assigned a quasar source with a luminosity value determined by our phenomeno- 
logical prescription. The location of the quasar from that point on is locked to the position of 
the dark matter particle closest to the center of the halo. The distribution of luminosities is then 
integrated from the higher value to the smallest up to the point where the average luminosity per 
unit volume reproduces the distribution fit given by Equation ([2]) at each redshift. For simplicity, 
we do not evolve the luminosity in each dark matter halo, which remains the same at initialization. 
As the luminosity function increases with decreasing redshift additional sources are spawned in the 
simulation, leading to an overall increase of their numbers. At z < 3 the flattening and subsequent 
decrease in the luminosity function is modeled by randomly removing point sources from the quasar 
list. In the right panel of Figure ([T]) we show the redshift evolution in the number of QSO sources 
in the volume. 

The location each point source is used to compute the 3D distribution of the tensor Eddington 



factor (§4.1) which is stored at each cosmic time step. The local value fij{x) at t" is an interpolated 
value between the local values of the two data dumps whose redshifts bound the instantaneous 
time (zi < t" < Z2). As we will discuss in the next sections, the cosmic evolution of the helium 
reionization is decoupled from the corresponding hydrogen one. Consequently, the placement of a 
He II ionizing source on a density peak with a precomputed hydrogen and therefore electron density 
can yield significant discrepancy between our calculation and a self-consistent one, particularly in 
close proximity to the source. 



3. Homogeneous Hydrogen Ionization 

As mentioned above, in the numerical and physical setup of this present work, we treat the 
ionization of H i and He i separately from He 11 . The metagalactic flux we use from Haardt &; 
Madau (2001) tabulates the contributions due to galaxies (GAL) and quasars (QSO) separately. 
In Enzo we have the option of running a simulation with GAL only, QSO only, or GAL+QSO. In a 
standard optically thin simulation the GAL+QSO UVB would be used. We have carried out such 
a simulation, hereafter called Simulation A, for comparison with the inhomogeneous reionization 
simulation. For the latter, we first run a hydro simulation using the GAL UVB to ionize Hi and 
He I . Then Hen reionization is accomplished by treating quasars as point sources as detailed 
in §4. Since our quasar population follows the same luminosity function and intrinsic spectrum 
as assumed by Haardt & Madau, we are able to compare the homogeneous and inhomogeneous 
simulations directly. We consider the effect of the e > 54.4 eV photon field as a perturbation 
on the previously ionized gas which affects only the He 11 -^ Hem chemistry while the Hi and 
He I ionization states remain unaffected. The treatment is by all measures an approximation that 
allows the problem to be solved as a single species ionization problem only and therefore it remains 
conceptually simple. In reality, the ionization of He 11 has two inter-dependent effects; it releases 
additionally one electron per He 11 ionization which in turn, when thermalized, raises the mean 
temperature of the IGM. These two effects combined would in principle shift the ionization balance 
of the H I /H II and He i /He 11 species. However, we show in this section that because by the time 
this takes place hydrogen and helium have already large ionization fractions the aforementioned 
effects are small. 

Starting with the rate equation for hydrogen and the cosmic mean density, we can safely 
ignore the collisional contributions and write the chemical balance in the proper frame of reference 
as follows: 

riHi = -3H{z)nHi - riniTi + nenHnai{T). (3) 

In Equation (l3|), rig = nnii + uhcI + '^nneii, Ti is the integrated Hi photoionization rate, 
and ai(T) is the radiative recombination coefficient in the Hii +e^Hi +7 reaction and is a 
function of temperature. In a cosmic medium, n^e = f^H, where f — 12 fo^ ^ mass of fraction 
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of ph = 0.75/9 and pHe = 0.25p respectively. The Heii number density can then be rewritten as 
nHeii = TT-He — ^Hei — ^ Hell I which in turn allows the electron density to be expressed as follows: 

rie = riHiXHll ^ Y2 J ^' 

Changing slightly the notation, we can rewrite Equation ([3]) according to Appendix (A) in terms 
of the ionization fractions X2 = riHii/nH, ipi = riHei/nHe, "02 = riHeii/nHe and Vs = riHeiii/nHe 
as: 

X2 = ri(l - X2) - X2nHai{T)[x2 + ^^(1 - V"! + V's)] => 

X2 ~ri(l-X2) -X2"'Hai(7')[X2H :^] (5) 

In the last equation, we assumed that almost all of the helium is highly ionized to the He 11 state 
ipi -^ 0. In ionization equilibrium {x2 = 0), after setting ^4 = (1 + ^3)/12 we can rewrite Equa- 
tion ([5]) as (1 + ^)j^-- = Ti/{AnHOci)- On one hand, the ionization of the He 11 would increase 
quantity A due to the increase in the Hem abundance. On the other hand, the increased tem- 
perature would decrease the recombination coefficient which can be approximated as ai oc T~^, 
for temperatures in the range of T ~ 10^ — 10^ K and /3 = 0.51. The approximation is based on 
expression fits by Bugress (1964). 

Therefore, we will investigate two extreme cases. Case (I) represents zero Hem abundance, 
■03 — > and would correspond to the unperturbed temperature T/. Case (II) represents a limit 
of almost full He 11 ionization, -(/^a — > 1, that would correspond to a new temperature T/j. In both 
cases, the hydrogen ionization rate is the same because our treatment of the individual QSO sources 
placement has a distribution that is statistically identical to the global average used for Hi and 
He I ionization. Therefore, we can form the ratio: 



{l + xi'/An)xi'{l-xi) ^ (Tn^pA^ , 

(l + ^yAj)xiil-x¥) ^Tj' An ^ ' 

Equation Q is solved in Appendix (B), for xi^ in the range of 1 — X2 = 10~^ — 10~^ and for 
Ajj = 1/6, Aj = 1/12, Tji/Tj = 1.5 — 2. The ratio of temperatures is based on estimates of the 
temperature increase due to the photoelectrons injected in the IGM from the ionized He 11 atoms 
(Haehnelt & Steinmetz 1998; Abel & Haehnelt 1999) and will be discussed in more detail in ^ 
The results in the Appendix figures show that when hydrogen is highly ionized then the shift in the 
ionization balance due to the He 11 ionization results in a decrease of the neutral hydrogen fraction 
primarily due to the decrease of the ionized hydrogen's recombination coefficient. 

The ~ 15 — 25% reduction depends on the temperature increase, where the largest increase 
(a factor of 2) produces the biggest shift. The fractional decrease in the neutral hydrogen frac- 
tion is smallest at low fractions of ionized hydrogen. For a typical neutral hydrogen fraction of 
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Xi — 10^^ at post-reionization redshifts, we can then estimate that our treatment (of not updating 
the hydrogen abundances) would overestimate the neutral hydrogen fraction by about 25% if the 
gas temperature increase is between 1.5-2 times the non-perturbed value. Consequently, we an- 
ticipate an overestimate by the same amount in the optical depth of hydrogen Lya radiation and 
an underestimate in the mean transmitted flux. The problem is similar to the one described in 
Zhang et al. (1997), Bryan et al. (1999) and Jena et al. (2005). An unmodified UVB by Haardt 
& Madau (1996), applied in those simulations, would not yield an agreement with the observed 
b-parameter of the Lya forest. An adjustment by a factor of 1.5-2.0 in the He ii photoheating rate 
was then necessary to match the observed results. Therefore, we caution against a strict interpre- 
tation of the resulting hydrogen Lya forest in the original calculation in which we suppressed the 
He II photoionization and photoheating processes altogether. 



4. Inhomogeneous Helium II Reionization 

We compute the inhomogeneous Heii reionization due to an evolving distribution of local QSO- 
type sources by solving for the time and spatial evolution of the ionizing radiation energy density 
Ey[T,t) as shown in following section. In §4.21 we describe our simplified chemistry model, and in 
r3]we present results. We defer a discussion of late photoheating to ^ 



4.1. Radiation Transfer Equation 

We consider simulation volumes of box length L much smaller than the horizon scale L <C Lh = 
c/H[z). Also, prior to bubble overlap, the time between emission and absorption of a random 
Hell ionizing photon will be much shorter than a Hubble time. In this limit, the cosmological 
radiative transfer equation reduces to the familiar one (Norman, Paschos & Abel 1998): 

— -KT H = riy- Xvh (7) 

c ot a 

where ly is the monochromatic specific intensity, r]u,Xu are emission and extinction coefficients, 
and v is the instantaneous, comoving frequency. In equation ^ the gradient is comoving, and 
hence we divide by the cosmological scale factor a to convert to proper distances. The zeroth and 
first angular moments of equation ([7]) yield 

— f + -V • F, = e, - ckyE, (8) 

at a 

and 

c^ ot a c 

where E^^Fy and F^ are the radiation energy, flux vector, and pressure tensor, respectively. All 
quantities are measured in the fluid (proper) frame, where e^ = 4:TTr]y is the emissivity, and k^ 
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is the absorption coefficient. The UV photons are scattered by the free electrons and for the 
range of electron number densities in the IGM (in the redshift interval we are interested in) we 

can estimate that Ag (the scattering mean free path) = > Lh- Therefore, we ignore the 

scattering coefficient. 

The radiation pressure tensor is coupled to the energy density Y'y — > Ej^ in the moment 
equations through the tensor Eddington factor, f^, = -g^. The latter guarantees the correct direction 
of the flux vector (Mihalas & Mihalas, 1984). It can be shown that the time derivative term in 
equation ^ is small compared to the rest if we integrate using a timestep long compared to the 
light crossing time of a computational cell, as we do. Therefore, dropping the time derivative in 
equation ^ and combining it with equation ([8]) we get 

^ = -V • [f V • {i,E,)] + e, - cKE, (10) 

where 

V.{i.E^)^^{rjE,). 

In Equation (jlOp . e^ is the spatially discrete monochromatic emissivity at the locations of the 
emitting sources and ky = nHeW^u is the local opacity, where the functional form of (J^; is 
given by Osterbrock (1989) and in Appendix C. Equation [10] can be solved for Eu(r,t) for a given 
source distribution provided the spatially-dependent Eddington tensor is known. Formally, fij is 
obtained from angular quadratures of the specific intensity (e.g., Hayes & Norman 2003). However, 
we wish to avoid solving the full angle- and frequency-dependent equation of radiative transfer. 
Instead we employ a geometric closure introduced by Gnedin & Abel (2001) in which we calculate 
the radiation pressure tensor assuming the medium is optically thin. In this limit 

p /N_ 1 Y^ Lk, {hi{r-Sk)){nj{r-Sk)) 

k 

where s^ are the positions of the ionizing sources, and hij are the direction vectors (basis) at the 
point r. Equation [11] describes pure radial streaming radiation from a collection of point sources 
in a transparent medium. Until He iii bubbles begin to overlap, this is an excellent approximation 
inside the He iii regions but a poor one outside. However, since there is very little ionizing radiation 
in the Heii regions, it makes little difference what one chooses for f. As discussed by Gnedin & 
Abel, the greatest error is when two bubbles begin to overlap and the two ionizing sources begin to 
"see one another." If one source is much more luminous that the other, this can lead to a ~ 10% 
error in the expansion rate of the smaller I- front. 

To solve Equation (jlOp we employ a finite volume method by rewriting the zeroth moment 
equation in a conservative form and integrating over a grid cell. The energy density, emissivity 
and opacity are zone centered quantities, therefore Jy E^dV = VgEy, Jy t^dV = Vge^, Jy kydV = 
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Vgk^, where Vg is the volume of a grid cell and Ey,ey, kj^ are now understood to be cell averages. 
Equation ([TO]) then becomes: 



dt aVg jce// 



BW 1 

"^ + -f7 </> ^u-d^cellsurf = (-v-ckyEy (12) 



dKy 0CCj 



Our time-implicit discretization scheme will be discussed in a follow up paper. Briefly, equation ()12p 
is discretized on a uniform cartesian mesh and integrated using backward Euler time differencing. 
Spatial discretization of the RHS of equation (fT2]) yields a 19-point stencil. The resulting sparse- 
banded system of linear equations is solved using the stabilized biconjugate gradient (BiCGstab) 
algorithm implemented in the MGMPI packagqj developed by the Laboratory for Computational 
Astrophysics (Bordner 2002). 

For this problem, we compute the radiative energy density at three frequency values above 
the ionization threshold value, Ei, E2, E^, corresponding to photon energies ei = 4, 62 = 2 x ei 
and 63 = 4 X ei in Rydberg units (e = ^^'^ ) respectively. The three points plus a fourth one 
at 64 = 32 Ryd with E4 = Es{y^)~°''' , are used to infer the interpolated profile of a 4th degree 
polynomial E{e)/Ei = Ylk=o^k{j)~'' between e = 4 — 16 Ryd. At e > 16 Ryd, we assume that 
E(e) = £^3(^)~"«, which follows from the reasonable assumption that at four times the ionization 
threshold energy, ~ 0.2 keV already in the soft X-ray energy band, there is little effect in the 
attenuation of the radiative energy due to opacity. The interpolation is necessary in order to be able 
to compute the ionization and heating rates which involve an integration in frequency space (photon 
energy) . The use of three frequencies and locally computing Ck bypasses the requirement more many 
frequency bins or the rewrite of the moment equations in frequency groups. Our motivation is that 
our calculation does not aim in computing the reprocessing of the radiation field spectrum, but 
rather in a reasonable and spatially inhomogeneous estimate of the photoionization/heating rates 
that will give rise to the 3D He 11 reionization process. 

Upon obtaining the solution -E1..3, E hereafter, we proceed to compute the local ionization 
and heating rates as described in ^4.2[ After that point, the coefficients c^ are no longer necessary 
and are discarded. The obvious advantage of this method is computational speed in the derivation 
of the photo-ionization/heating rates through an analytical formula. The disadvantage is that the 
accuracy is as only good as the cubic interpolation scheme. However, we note that our interpolation 
scheme does a reasonably good job of describing the processed quasar spectrum obtained with the 
full multifrequency calculation of Abel &: Haehnelt (1999). In fact, our choice of frequency points 
was strongly guided by the inset spectra in their Figure 1. 

The distribution of local sources determine the point source emissivity (source function) and the 
3D distribution of the Eddington factor. Because in our numerical setup all sources emit radiation 



^Ica.ucsd.edu/portal/software/mgmpi 
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with identical spectral slope Uq, the calculation of the Eddington tensor via equation (jlip yields a 
quantity that does not depend on frequency {fu = p^). The pre-processing involves keeping track 
of the sources position and luminosity during the calculation and updating the Eddington factor 
functional form and emissivity source function at the beginning of each data dump calculation. 
The Eddington factor is initialized at flource = l/SiJ*-' within each grid cell containing a source. 
Between time steps, assumed to be the redshift interval between the data dumps {5z = 0.1), we 
solve for the E{r,tn) which we use to update the Heii photo-ionization rate F™ described in §4.21 
The photo-ionization rates are then used to compute the next time-step for solving the transfer 
equation. Because the latter is solved implicitly, the solution is not sensitive to any particular 
choice of the time-step, beyond obvious concerns of numerical convergence. 



4.2. Chemistry Implementation 

Our simple single species chemistry determines the time step of the evolution. We update the 
ionization fraction ?/;3(x, t) = — ^^^^ through the rate equation: 

nHelll = -3H{z)nHeiii + nHeII^2 - nenHeiiia2(T) (13) 

In Equation (J13p . r2 = T™ -|- T2 is the sum of the radiative and collisional ionization rates. The 
collisional ionization is due to collisions of the He 11 ion with electrons. Hen + e — > Heni + 2e, 
and therefore is proportional to the electron density ng. The collisional ionization rate per electron 
Eg /rie is a function of the gas temperature. An analytic fit to the temperature dependence is 
provided in Appendix C where we also provide the functional form for the recombination coefficient 

(both quantities are measured in cm?s~^. The photoionization rate per baryon is then given 

00 

by the equation (Osterbrock 1989) F™ = J cj^a^ dv. Because we find a parametric fit of 

the energy density in frequency space the photoionization can be directly computed as follows. 
^2"'^ = Hell ^[-^3 0^3 + ^'^Ylik=o 'i+^{^ ~ 4"^'^'*'^^)], where we approximate the ionization cross 
section with the power law 0"^ = o"224 (i)~^- 

To determine the time-step (5t" for the advancement of the radiation solution between f" and 
4"-+^, we follow the procedure below. We collect the photoionization time-scales, Tchem = (E™*^)"^, 
from grid cells that lie in the vicinity of the ionization front and then calculate 5t through 5t = 
max {Tchem, Tic), where Tic is the cell light crossing time-scale. Tic = 5xc~^. Cells close the I-front 
interface are "captured" by the criterion ^""'^ (^") ^ 0.1. Comparison between the light-crossing 
and chemical time scales is necessary in the initial moments of the evolution, because in proximity 
to the source the I-front propagates close to the speed of the light. Evolving the energy density with 
the light-crossing time-scale constrains the I-front expansion to subluminal speeds. In addition, the 
use of a chemical time-scale close to the source would yield very small time-steps, due to the large 
number of photons, that would in turn slow down the overall calculation. 

Before we proceed, we need to make the following very important clarification. In this 
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work, we do not consider contributions to the radiative energy density by Hell recombinations 
despite including them in the update of the Heii abundance. The significance of the contributing 
Hell recombinations to the energy density are important when a consistent ionization calculation 
of all species is performed. In that case, radiative recombinations to energies e < 54.4 eV would 
contribute to the ionization balance of hydrogen and helium. However, these species are already 
at high ionization fractions and therefore such contributions are not expected to have a significant 
effect in our scheme. In a consistent calculation, the photon fiux from the central source creates a 
concentric set of Stromgen regions where the ionized hydrogen and helium II extend ahead of the 
corresponding helium HI volume. In such a scheme, recombinations at the Hem ionization front 
yield photons with energies e < 54.4 eV which can freely propagate forward through the already 
ionized Hei and Hi and assist in the advancement of the Heii and Hii fronts respectively. In our 
setup, there are no He ii and H ii Stromgen regions, only He iii ones. Therefore, the photon fiux 
from recombinations on the He III I-fronts are ignored and we only consider the chemical effects 
resulting from the attenuation and percolation of the individual He ii ionizing flux. 

However, recombinations are included in the abundance update primarily because it can be 
an important effect for the diffuse helium in proximity to a local source that shut off. The lack of 
direct photons could lead to a rapid recombination of the He iii bubble, if no additional radiative 
fiux reaches that region from another QSO source. Equation (J13p can then be rewritten as follows. 



V-s = ^^ ~ ^^^ - ^ (14) 

In Equation ([H]), Tg"" = (r2)^^ is total ionization time scale and rl^*^ = {a2{T)ne)^^ is the local 
He III recombination time scale and r2 is the sum of all ionization processes that lead to the forward 
reaction Hell — > Hem and we have assumed V'l ^ 0. These processes in our scheme involve the 
combination of the direct photoionizations from the QSO sources F™'^ and collisional ionization 
f.k'^°\T). We integrate equation 1141 using backward Euler time differencing where all source 
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terms are computed at the advance time t 

^n+1 ^ i^^+dt/[T^ ) 

In Equation (jlSp . only the photoionization rate r™*^ is available at the advanced time. There- 
fore, we are forced to initialize the abundance update by computing (:;^)"+^ '^ (r'-'"')"+i+ng(fc^°'(T))" 
and substitute (:p4c)"''^^ -^ (:;:4c)"'. However, upon obtaining the Hem number density at t""^^, 
iT'Heiii = V's^ f^He we Can update the electron density at t""''"'^ and insert it back in Equation (|15p . 
We therefore can improve upon the original estimate by iterating Equation (J15p until — ^ < 0.1. 
Unfortunately, our methodology of updating the temperature, as described in ^ is crude and is 
not used in the iterative scheme. The local helium number density at any time is computed from 
the gas density as riHe = £^ — T6m~' "^^^ electron density is given by the charge conservation 
equation rig = riHjj + riHeii + '^nneiii- Our calculation assumes no change in the ionized hy- 
drogen density between the value in the original simulation (°), and the post-processed value (^). 
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Therefore, n\-nl = {n\^^jj - n%^^^) + 2{n]j^jjj - n^^jjj) = {n\j^jjj - njj^jjj). which allows the 
calculation of the local electron density at any time from iie = n° + {nHelll — ''T'Helll)- 



4.3. Results 

In Figure ([2]), we show volumetric renderings of riHeiiii^, z) at two redshift instances. The 3D 
visualization shows the expanding ionized bubbles filling up the cosmic volume due to the combined 
effect of the radiative energy transport and sources being turned on at different parts of the volume 
at later redshifts. Individual bubbles of Hem may stagnate as they reach their Stromgren radii 
due to recombinations. However, overall the volume filling factor (VFF) of Hem increases as 
more quasars are placed in the computational volume and the percolation between the I-fronts 
increases the mean free path of the ionizing photon flux. In left panel of Figure (l3|), we show 
a slice through the cosmic volume at z = 2.6 of the Heii , Hem density distributions. Ionized 
regions have percolated through the cosmic medium to "open up" the IGM to > 54.4 eV radiation, 
effectively completing Heii reionization by such redshifts. In the right panel of Figure ([3]), we show 
the redshift evolution of the VFF, as measured by the fraction of the grid cells with ionized helium 
at Hem abundance of ^3 > 10~^. As the ionized regions begin to merge, assisted by the increase 
in the QSO number density, the VFF(z) rapidly increases, leading to a value of > 68% at z < 2.8. 
The redshift of significant merging, which we define as VFF ~ 0.90 is achieved by z ~ 2.5 where we 
point to a statistical global He 11 reionization. This redshift value compares well with the observed 
determination by Kriss et al. (2001). The solid line in the VFF evolution figure is a spline fit 
through the computed data {5z = 0.1). For reference, we include the derived values every 5z = 0.5, 
along with the error estimates based on the location uncertainty of the I- front. The uncertainty is 
simply due to the fact that computing the radiative energy density in the zone centers yields no 
information on the profile of the field across the grid cell. Therefore, we assigned an error estimate 
in the ionized volume fraction equal {^Y /V ^ where 5x is the grid resolution and V is the volume 
of the computational box. The evolution of the VFF shows a rapid increase in He iii at redshifts 
z < 4, following an earlier epoch of apparent stagnation. 

An alternative way to illustrate He 11 reionization is to plot the redshift evolution of the volume 
averaged abundance fraction. In the left panel of Figure ([!]), we show that the mean mass fraction 
in Hell puenl P^ drops significantly at 2: < 4. For reference, we also show the mean fraction in 
Hi phiI P-, undergoes a steep drop at z ~ 6.5 and continues to decrease under a smooth redshift 
profile, the properties of which are discussed in Paschos &: Norman (2005). One notable difference 
in examining the two evolution profiles emerges between the slopes of the two curves. Beyond 
the fundamental differences in the ionization calculation of the two species, the result shown in 
Figure ([!]) shows that the He 11 reionization epoch is much more extended than that of H i . A 
rapid drop in the Hi fraction occurs between z = 7 — 6.5 (Razoumov et al. 2002; Paschos & 
Norman 2005). In that redshift range the mean Hi mass fraction drops by 4 dex. By visual 
inspection of the mean fraction in He 11 from Figure ^ we can determine that it drops by ~ 2 dex 
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between z = 3.5 — 2.5. When the difference in the redshift interval is converted to cosmic time 
it yields that a similar drop in the mass fraction takes approximately eight times longer to occur 
in the case of He ii when compared to H i . This is consistent with the conclusion reached in ^2.2^ 
where the latency in the Heii reionization is attributed to higher recombination time scales and 
less available ionizing photons per baryon for a stellar dominating UV background. 

In the right panel of Figure (JH), we plot the Heii mass fraction versus the local overdensity 
by logarithmically binning the latter quantity and computing the mean and median He ii from the 
cells with gas density within the bin. Such a graph is intended to show the trend between the two 
quantities. For reference, we plot the trend between the two quantities from a standard cosmological 
simulation, where all ionizations are computed self-consistently due a uniform UV background in 
the optically thin approximation. 

The simulation is terminated at z ~ 2.5 at VFF ~ 0.9. The reason for suspending the 
calculation is that when the cosmic volume is effectively transparent to the ionizing radiation, 
the assumptions underlying equation ^ no longer apply. Specifically, free streaming photons 
may cross the volume unimpeded by absorption and their mean free path can become comparable 
to the size of the horizon. The latter effect requires keeping the cosmological dilution term in 
Equation (112p . which was ignored. From the numerical perspective, solving a parabolic equation 
when the conditions call for a hyperbolic one, can create local superluminal speeds of the ionizing 
front if the chemical time scale becomes smaller than the cell light crossing time. Concluding the 
calculation at the end of the opaque phase of He ii still addresses the main question investigated in 
this work; what is the primary mechanism that leads to the He ii reionization. We conclude, that a 
rising population of QSO sources, assisted by the gas dilution due to cosmic expansion, dumpiness 
due to structure formation and the pre-ionization of neutral hydrogen which allows for the almost 
exclusive usage of the He ii ionizing radiation are a set of physical conditions that reproduce the 
epoch of Hell reionization by redshifts z < 2.5. 



5. Late Heating due to He ii Reionization 

5.1. Physical Considerations 

Photoionization of Heii at an epoch later than that of Hi releases an additional one photo- 
electron per ionization. We can estimate the mean energy of such electrons due to absorption of 
photons with energies > 54.4 eV= 4 Ryd by He ii and compare it to the value obtained in the case of 
> 13.6 eV= 1 Ryd absorption by Hi atoms. For simplicity, we will ignore the geometric attenuation 
and optical depth effects in the radiative flux due to local sources and gas opacity and assume that 
the local mean intensity of the radiation is the same as the emitted, with a spectrum J^ = jQi2e~°'i. 
As before, in this notation e = ^^'^ . The photo-heating rate, in ergs/s, due to electrons ejected 
from Hell atoms, can then be computed from Gueii = ^912/4 (4vrJe/e)(e — 4)aede. Substituting 

for the power law of the mean intensity, for a^ = '^He//'^^^"^ ^^ Set Gueii = ^'^"(q '^l^v'^"'!^) • 
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In a similar manner, we can derive Ghi = (q. +2)(a +3) ■ ^^ ^^^ sources responsible for the Heii and 
H I ionization are the same then the J912 amplitude and the Oq spectral slope are identical in both 
equations which allows for the derivation of the following ratio: 

Ghi o-hi 

Substituting for o-|^g// = -f^ we get ^h„ii = 4 09 The total photoheating rate per unit 
volume is proportional to the number density of the absorbers which yields the ratio between 
He II and H i photoheating rates to be Quell I Ghi = "'^"^^ 4""'? . The ratio of number densities is 
estimated in §2.21 to be nHeii/nui = ^S, where S is the softness parameter. The ratio of the 
photoheating rates then becomes GhbIi/Ghi = ^5*4""" ~ A4"9+i4-Qg = 5/3. The latter factor 
is indicative of the degree of temperature increase due to the thermalization of the photoelectrons 
ejected by He 11 ionizations when compared to the temperature inferred by HI ionization alone. 

In deriving the above value, we made the assumption that hydrogen and singly ionized helium 
are photoionized simultaneously by the same ultraviolet field. The effects of distinct ionization 
epochs can however be modeled in the above relation if we adopt the scenario of hydrogen ionization 
by a soft radiation background {a^"^' ~ 5) and of He 11 ionization by a hard one {a^^' ~ 1.5). At a 
point in time when the populations of galaxies and QSOs have the same energy output at the LyC 

(2) 

limit (z ~ 4), we can derive an estimate of the photoheating rates ratio to be q'^" = f^^-fi) — - — -j 

for Oq = 1.5 and for Og =5. In conclusion, the thermalization of the photoelectron in the 
Hejj + 7 — > Hejjj + e~ reaction can be an important determinant of the intergalactic medium 
temperature. 

Photoionization models based on hydrogen ionization alone predict a temperature of the inter- 
galactic medium of Tjcm ~ 1-2 x 10^ K. However, Lya forest observations yield lines with median 
broadening widths of 6 = 26 — 36 km/s (Carswell et al. 1987, 1989; Zhang et al. 1997; Dave et al. 
1997) in the intermediate redshift range of z=2-4. Because the thermal and differential Hubble flow 
components in the total HI line broadening are of the same order of magnitude (Zhang et al., 1998; 
Aguirre, 2002), one can infer a thermal width range in the HI Lya forest between bth = 13 — 18 
km/s. The effects of the peculiar velocity can be ignored if the the focus is at densities close to the 
cosmic mean. This value range in bth would require the temperature in the intergalactic medium to 
be TjGM = 20, 000 — 40, 000 K, a factor of ~ 1.7 — 3.3 above the temperature inferred by hydrogen 
photoionization alone. As seen above, such an increase in temperature can be reproduced by the 
Hell ionization photoheating due to a hard ultraviolet spectrum. 
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5.2. Perturbed Gas Energy Equation 

The post-processing of the simulation data dumps involves an update of the local gas temper- 
ature by solving for the perturbed temperature in the thermal energy equation. In the simplest 
approximation, we will assume that the change in the net rate of the cell thermal energy density 
is due to the balance between the heating by the ejected photoelectrons from Heii ions and the 
Hem recombination cooling, 7:—ij;^^-st = G — A, where 7 = 3 is the gas adiabatic index, p the 
local gas density and /i the mean molecular weight. All quantities have local proper values. In 
addition, in the notation followed here G = nneiiGr is the total photoheating rate, measured in 
ergs/ s/ cm?, which is the radiative rate Gr = G mentioned above. The cooling rate. A, is a function 
of the local gas temperature and density and is equal to A = nenHeii[L{T). L denotes the cooling 
function due to Hem recombinations and is given by L{T) = 3.48 x 10-'^^T^/^T^^-'^{1 +T^-'^)-^ in 
units ergs cm^/s (Cen 1992), where T„ = T/IO". 

Recombination Hem cooling is only one out of several processes that cool the cosmic gas, 
the list of which is described in detail in Anninos et al. (1997). The cooling rate coefficients, 
in parametric fits that depend on the local gas temperature, due to excitation, ionization and 
recombination of the primary chemical species along with bremsstrahlung and Compton cooling 
are fully incorporated into the ENZO code. We expect that Hem cooling would dominate the 
cooling processes inside the He iii bubbles when compared to the corresponding He 11 recombination 
cooling, because the He 11 abundance is significantly reduced there to fractions ^2 — 10^^ — 10^^. 
In addition, the Hem recombination cooling rate, along with Hii recombination, dominate at the 
low temperatures found in underdense cosmic regions T < lO'* K over all other types. In a scheme 
where hydrogen is already almost completely ionized at 5 < 1 {x2 ~ 1); but helium is predominately 
only singly ionized ('(/'2 ~ 1) the thermal balance would be controlled by the heating and cooling of 
the latter species' ionization. 

However, our ability to accurately post-process the temperature field in our scheme is limited 
by the fact that excitation and collisional ionization cooling from processes involving H i , He i and 
Hell species dominate the cooling curves at temperatures T > 10^ K, while exhibiting very steep 
profiles in the temperature range T ~ 10^ — 10^ K. An increase in the gas temperature by the 
photoelectrons ejected in the He 11 radiative ionization would cause an increase of the cooling co- 
efficients. Combined with the increase in the electron number density this shift in the thermal 
balance could significantly increase the cooling rate even though the fractional abundances in Hi , 
He I and He 11 are small. The end result is that, by only including Hem recombination cooling in 
recomputing the gas temperature, we may overestimate the value of the adjusted temperature. This 
upper limit in the temperature estimate implies that the recombination time-scale in Equation (J15p 
is also an upper estimate and that the overall propagation of the cumulative He iii ionization front 
is faster than in the case of a self-consistent calculation. However, such calculation would require 
computing the ionization and thermal balance of all species self-consistently and is beyond the 
scope of this paper. A numerical scheme is under development that will allow us to do this in the 
near future (Reynolds et al., in prep). Therefore, we conclude that the epoch of He 11 reionization 
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may be placed at a later redshift than the one we calculated here {zreion ~ 2.5) even though we 
note that such adjustments could be reversed or may not be necessary if the diffuse recombina- 
tion radiation is added in the calculation. Such radiation would further ionize HI and Hel species 
suppressing their contributions to the cooling curve. That may explain why, even under all of the 
assumptions and approximations that we allowed and followed in this work, the end result of our 
predicted redshift evolution of the He ii opacity correlates well with observations of the He ii Lya 
forest as we shall show in ^ 

We proceed with further detailing our temperature update method. If T^"' and T^^' denote the 
cell temperatures before and after the presence of Heii ejected photoelectrons, then we approximate 
the change in the thermal energy equation as follows: 

J^.^+'^-V"''""'(GW-AW) (17) 

ot ot kp 

In an explicit, time-discretized form, where the original temperature at time i" is obtained by 
interpolation between the logarithm of the local temperature in the two data dumps with redshifts 
that bound the time evolution {zi < t" < Z2), Equation (fT7|) becomes: 



r«i - T^i) = Ti% - tW + {6tUj - 1) 

li^r^HnthGr'^ - nf^n^lL-) (18) 

In Equation (jlSp . time-centered quantities are computed as X"'~''2 = 0.5(X"+^ + X") and 
only for variables that we know their value at the forward time i"""*"^. Therefore, L is computed 
at time t" because it is a function of temperature which is unknown at t"+^. The update in 
temperature occurs after all principal quantities of i?"'^^, nfjeii-, "n-Hem and He are computed at 
^n+i Y)Y advancing the local solutions at the implicit time-step of the radiation field discussed in 
§4.21 The quantity ^^^^^ = ^ rij is equal to total number density of the cosmic species plus electrons 
and is computed as follows: ^ = [^:^ + n^'^ ~ [0.9 • 10"^J^6/i^(l + zf + n^'^ , where the 
expression for rig was provided in ^ 

In an ideal calculation, if a local grid cell is outside the Hem bubble at f^ then T„ = Tn 
and n^g — > 0. If no direct ionizing radiation reaches that grid cell by t*^"*"^ then n^ -^ and 
Equation ([TH]) would predict T\^J^^ = T)^J^^ . In a cell where thermal equilibrium was reached during 
the original calculation between t^~^^ and t" then during post-processing T^J.i = r„ is achieved 
only when the heating and cooling terms balance out. 

A point of concern is that there was no physical reason to suppress the collisional ionization of 
Hell in the original simulation. Even if there are no ionizing sources capable of radiatively ionizing 
Hell , a local temperature of T > 4.64 • 10^ K, found in overdense regions, may be enough to 
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collisionally eject the 54.4 eV bound electron in the Heii atom. Consequently, the temperature 
evolution T^°'(z) includes such an effect and can skew the post-processed evolution T^^'{z). This is 
evident if we assume that due to collisional ionization in the pre-processed data nHem 7^ and at 
r, Xl^^ = Xi°^ Equation ^ would then predict TJ^% = T^J^ + {6t)nW{n^°ljj,n^°l^^^,ni°\T(°)) 

where W ^ (^ - l)l(!^)"+i x (-1) n:J^4+J^^^^L"(r(°)) 

The last equation would unnecessarily recompute the temperature in regions which lack ra- 
diative input but have significant collisional rates in Heii <-> Hem and therefore, nHeju,o 7^ 0. 
To account for such discrepancy in the temperature evolution, each updated local temperature 
is adjusted as T^_li = T^J.i + ('^*)nW^('^^e//' '^He 5''^e i^ ■*)• Although in doing so we improve 
upon the temperature evolution, the collisional effect can never fully be readjusted because of 
interdependency between all of the physical quantities. 



5.3. Results 

In Figure [5l we show the results of our temperature calculation. In the left panel, we plot the 
evolution of the mean temperature in overdensities log{5) =0 — 1 (solid curve). Overplotted is the 
evolution in the original calculation (dashed) to showcase that at z < 3.5 the two profiles are devi- 
ate. This demonstrates the effect of late Heii photoheating due to the rising population of QSO's 
in the cosmic volume. The overdensity interval was chosen in order to show that He ii reionization 
is a cosmic event that primarily affects the diffuse and mildly overdense IGM where the the tem- 
perature can get increased by about a factor of 2 at the end of the calculation. One should consider 
two effects that support this conclusion. In regions of significant overdensity and therefore dark 
matter potential, gas is heated by the gravitationally controlled free fall compression to temper- 
atures T > 10^ K. Therefore, the effects of photoheating due to the photoelectrons ejected by 
radiative ionizations of Heii only result in an insignificant fractional change. In addition, the large 
electron density and gas density can be a source of large optical depths and radiation trapping 
due to increased recombination time scales for a modest increase in the temperature. As a result 
the radiative energy density can significantly drop in such regions which furthermore reduces the 
efficiency of He ii ionization. 

The optical depth effects are demonstrated in the right panel of Figure ([5]), where we show 
(z=2.5) the median temperature vs. density (solid line) in the scatter plot between the two quan- 
tities on the simulation grid. We obtain the curve by binning the gas overdensity in logarithmic 
intervals and computing the median gas temperature within each bin. For reference, the relation 
in the original calculation is also shown (dashed lines). On the left panel of Figure Q, we plot the 
ratio between the optically thick optically thin calculations. The dashed lines shows the effect of 
our postprocessing on the gas temperature. Heii photoionization contributes primarily to the gas 
temperature at the mean and low gas densities. The effect is diminished at higher overdensities 
where collisional ionization dominates. In addition to comparing to the original calculation we 
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also compare to a simulation where the chemical species abundunces are self-consistently computed 
during the simulation due to the same homogeneous UV background and in the optically thin limit. 
The dot-dashed curves on the right panel of Figure ^ and the left panel of Figure ^ trace the tem- 
perature density relation in that case and the ratio to our optically thick calculation respectively. 
We find that our postprocessed temperature is higher by a factor of ~ 1.7 at the cosmic density 
level when compared to the optically thin calculation with a homogeneous UV background. This 
is consistent with the analytical approximation between an optical thick and optical calculation 
derived in Abel & Haenhelt (1999). Finally, we plot on the right panel of Figure ^ the relation 
between the slope of the equation of state versus the gas overdensity for all three simulations. 

In conclusion, the increase in temperature is a manifestation of the additional heating that 
results from photoelectrons ejected by Heii ionizations in the redshift interval that corresponds 
to the rise in the number density of sources emitting hard radiation. The physical implication is 
that the late increase in temperature may well be the reason why the galaxy luminosity function 
decreases at z < 4. A raise in the mean temperature due to the cosmic evolution of QSO sources 
would increase the Jeans mass threshold by a factor of 2.2-5 if we adapt an average increase in 

3 

the temperature between 1.7-3 {M jeans oc T2). The latter would in turn suppress the further 
formation of dwarf galaxies in the cosmic volume. However, we find that, according to the right 
panel of Figure ([5]), the fractional increase in local gas temperature is on average not that large at 
higher overdensities. If we exclude cosmic neighborhoods that are in close proximity to the local 
UV sources, where the temperature increase is significantly greater due to the high radiative energy 
density, then in collapsed structures at 6 > 100 which are increasingly self-shielding, photoheating 
is not as effective as shock heating due to gravitational collapse. Nonetheless, a firm conclusion 
on that effect is not possible in this work, due to the very coarse grid resolution that does not 
adequately resolve the aforementioned structures. 



6. Signatures of He II Reionization 

6.1. Synthetic Flux Spectra 

The updated He 11 density and gas temperature can be used to study the effects of this in- 
homogeneous reionization scheme on the transmission of the intergalactic medium transmission in 
the Hell Lya restframe wavelength. The objective is to compare with the He 11 transmissivity ob- 
tained from analyzing the currently available observed lines of sight. We synthesize spectra of the 
Hell Lya absorption at the rest wavelength of 304 A along 300 random lines of sight (LOS). The 
number of LOS was chosen in order to yield a less than 1% fluctuation to the mean transmitted 
flux by the end of the calculation. 

The synthesis method is described in detail in Zhang et al. (1997). In addition to He 11 , we 
also compute the transmitted flux of the corresponding Hi forest at the rest wavelength of Lya 
1216 A. Each velocity pixel registers the local absorption in He 11 and Hi and therefore maps onto 
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the same grid location for the two redshifted wavelengths. The redshift interval of the spectrum 
output is set at Az = 0.2. In addition to this interval being the redshift output interval of the 
hydrodynamic simulation it is a long redshift path that minimizes the sightline to sightline variance 
by forcing the transmission path to cross and wrap through the volume boundaries along the same 
directional vector about ^^ ^ times. /\zcuhe = L H{z)/c is the linear size of the cube in redshift 
units where H(z) = 100/i((l + z)^Qm + ^a)^ is the Hubble constant at redshift z and L=100 Mpc 
comoving. At z = 2.5 and for h=0.71 the Hubble constant is H{z = 2.5) = 249.1 Mpc/km/s which 
yields Azcube = 0.083 and yr — ^— = 2.4, the number of wraps through the computational volume. 

The output from the hydrodynamic simulation is saved every dzdump = 0.2 and each optical 
depth integration is computed between two data dumps at zi and Z2- The output spectrum is 
therefore centered at 0.5 x (zi + 2:2), where Z2 = zi — dzdump- The input fields are species density 
(Hi and He 11 ), gas temperature and the three peculiar velocity components, all assumed to be 
frozen in the comoving frame of reference. However, we allow proper evolution along the redshift 
path of the sightline for the densities (oc (1 + 2)^) and velocities (oc (1 + z)~^). The spectra are 
computed along lines of sight that sample the computational volume continuously for z < 6.1. We 
do this by restarting the calculation in the new redshift interval from the mesh location where 
the previous calculation stopped and continue along the same directional vector. The initial point 
of origin is randomly selected. After the completion of each Az = 0.2 segment we can paste 
all segments together to obtain a continuous transmission line of sight. The latter can then be 
resampled at intervals centered to suit our analysis. 

At redshifts where absorption features can be recognized as blends of individual lines, a de- 
blending algorithm, described in Zhang et al. (1997; 1998) based on fitting lorentzian profiles below 
a transmission cutoff, can identify such lines and compute properties such as the column density, 
broadening width (b-parameter) and equivalent width. We must note, that as the opacity of the 
IGM increases with redshift and the transmission spectrum is dominated by dark regions, the al- 
gorithm fails in identifying the lines. Typically, good results are obtained just past the reionization 
opacity tail which will be the focus in this analysis. 

In Figure d?]) we show an example sightline of H i (top left panel) and He 11 (bottom left panel) 
Lya transmission at z=2.5 and compare with the results from an Enzo simulation that computes 
abundances self-consistently in the optical thin limit (top and bottom right panels) using a UVB 
from a mix of evolving quasar and galaxy populations (Haardt &: Madau 2001). Both transmis- 
sion lines are casted from the same initial point and along the same directional vector. We use 
the updated temperature and the unperturbed neutral hydrogen abundance for calculating the 
Hi transmission spectra. Although there are few visible differences at first glance, the Hi spectrum 
in the postprocessed temperature case has a continuum flux level below that of the optically thin, 
self-consistent case. Since the mean transmitted flux is sensitive to the number of pixels close to 
the continuum at low redshifts the line of sight average transmitted flux in the post processed 
calculation is below the value obtained from the self consistent calculation by about ~ 8%. That 
was expected in our discussion in Section ([3]). 
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In the absorption spectra along the randomly casted sightlines there is an apparent lack of 
signatures marking the presence or not of the local ionizing sources. The spectra in the examples 
shown in Figure ([7]) look in general similar, although in some isolated regions the local amplitude 
clearly differs from a uniform UV background. This result is expected because at the time of 
reionization completion, the mean absorption in a random sightline spanning 2.4 times the length 
of the simulated cube is not sensitive to any local effects by the ionizing sources but rather depends 
on the value of the ionizing radiation at the level of the mean density and the density valleys. 

We show on the left panels of Figure ([8]) a sightline with a path that takes the sightline about 
one mesh resolution element, ~ 0.2 Mpc comoving, away from a UV source. On the top left panel, 
we show the Lya He ii transmission versus observed wavelength. Below, we plot the instantaneous 
comoving distance of each velocity pixel to the closest UV source versus the comoving path length of 
the sightline. The last curve is a joint ensample of parabolic curves each having a minimum, marked 
with crosses, that corresponds to the position of nearest distance to the ionizing source. Each time 
another UV source becomes closest to the trajectory, the curve is marked by the beginning of an 
another parabola. We note that at the location of closest proximity for that sightline, at ~ 120 
Mpc comoving along the path, the corresponding transmission is almost 1. Since the sources are 
placed at the high density dark matter peaks, which they subsequently follow during the dynamical 
evolution of the dark matter distribution, the absorption at that location is due to the UV source 
proximity. However, at a part of the spectrum about 1042 A, the high transmission there is due to a 
highly ionized underdense region since the closest point source is more that 20 Mpc comoving away. 
Therefore, from the information in the spectra alone, we will not be able to distinguish between 
transmission in a region close an ionizing source or transmission from an underdense region. 

We will call impact parameter each nearest distance to an ionizing source encountered by 
the sightline trajectory. On the right panel of Figure ([8|), we plot the sightline distribution of 
impact parameters in bins of constant size 1 Mpc. The upper horizontal axis converts comoving 
Mpc distances to proper separation velocities. The distribution is negatively skewed with a peak 
frequency (mode) located at ~ 9 Mpc (~ 550 km/s). Due to negative skewness, the mean is shifted 
to the right at a value of ~ 14 Mpc (~ 750 km/s). The negative skewness is most likely due to 
the fact that the distribution of the UV point sources on the grid is not isotropic but are clustered 
according to the clustering properties of the host dark matter halos. From the distribution in the 
right panel of Figure ([8]) we can compute a typical range of impact parameters in the sightlines of 
~ 2 — 22 Mpc comoving contained within the values of the distribution at 1/e the peak value. We 
are interested to determine whether the transmission properties of our sightline sample is in any 
way biased by the proximity to the ionizing sources or if we are mostly sensitive to the values of 
the ionizing flux at the mean level. 

Detailed analysis of the proximity effect around each source is beyond the scope of this paper. 
However, we can impose an upper bound limit based on the highest luminosity quasar in the 
computational volume at z=2.5, L™*^^ = 2 x 10'^^ ergs/s where L4 designates continuum ionizing 
luminosity above 4 Ryd. An estimate of a proximity size distance can be obtained from the 
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volume averaged energy density above the Heii ionization threshold, -E4. In our calculation we 
have computed E4 = 7.7 x 10~^^ ergs/cm~^ which corresponds to a volume averaged He 11 ionizing 
flux equal to -Fli = cE^i = 2.3 x 10~^ ergs/s/cm~^. The condition that the geometrical attenuation 
of the central emission equals the total flux through the surface of a sphere centered at the source 
yields an upper bound for this proximity distance of (f^'^Y ^ {L'J^"''^ / ^1:^^)2 ~ 0.86 Mpc proper 
or dfjeii ^ 3 Mpc comoving. When compared to the comoving length of the sightlines (~ 270 
Mpc comoving) we expect that these regions do not have a large effect on the mean He 11 Lya 
transmission value. 

On the left panel of Figure Q, we show a quantitative estimate of that effect. There, we 
scatter plot the average impact parameter of the sightline versus the mean transmission in each 
sightline (cross symbols). The degree of their linear correlation, measured by the Pearson correlation 
coefficient, is rather small, r = —0.24. The figure does show that there is non-negligible negative 
correlation between mean transmission and the average proximity to the distribution of ionizing 
sources for each sightline. A power law fit to the data in the form Fios oc d~^ where Fios is the 
sightline specific mean fiux and dip is the average impact parameter yields s = —0.16 it 0.04 and is 
plotted over the individual points. The curve lies within the two standard deviations levels of the 
mean transmission. Therefore, we conclude that the proximity effects have no statistical weight on 
the sightline transmissions. That is consistent with our estimate of the upper limit on the proximity 
distance (~ 3 Mpc comoving) based on matching the attenuation of the point luminosity to the 
volume averaged ionizing flux. Since the typical range of the impact parameters is between 2-22 
Mpc comoving our sightlines are on average too far away from the point sources. This conclusion 
is also supported by the right panel of Figure ([9]) where we show the scatter plot between the 
individual impact parameter values and the local He 11 Lya absorbed flux at that distance from 
the point source (points). The data show a tail of low absorption at small impact parameters at 
distances < 4.5 Mpc comoving. The distribution of absorption above > 4.5 Mpc is consistent with 
the mean level of absorption in the sightline sample, shown as a horizontal blue solid line along with 
upper and lower blue dashed lines corresponding to the 2a level. To get a clear view of the trend 
we bin the horizontal axis into constant logarithmic bins and compute the median (red histogram) 
and mean (blue histogram) per bin. A Fermi-Dirac function in the form /(x) = tdo-x) is then 

l+exp{- - 

fitted to the median histogram with do = 3.1 Mpc (point of half maximum) and /i = 0.63 Mpc (skin 
width) for impact parameters up to 10 Mpc. The point of half maximum is consistent with the 
previous estimate of the proximity effect distance based on the luminosity and mean intensity of the 
UV field. The figure clearly shows that the effect of source proximity on the local absorbed fiux are 
relevant at comoving distances of < 7 Mpc where the local absorption matches the mean He 11 Lya 
absorption in the IGM and becomes significant at impact parameters < 4.5 Mpc comoving where 
the local absorption matches the lower 2a level of the mean absorption. 
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6.2. Optical Depth Evolution 

The straightforward average of the flux, F = jv x-nios ^l^°'^'^i ^^-^«i' ^o^^ ^^ pixels and lines of 
sight per redshift interval, is a measure of the opacity of the cosmic volume due to the absorption by 
the particular chemical species (Heii and Hi in this case). In the notation followed, Npx = 30, 000 
is the number of pixels per redshift interval, which yields a redshift pixel resolution of Rz = 5z = 
0.2/30,000 ~ 6.6 X 10~^, and nlos = 300 is the number of random lines of sight. The redshift 
resolution is fixed in our calculation, which results in a redshift dependent spectral resolution of 
i?A = AA ~ ^JT ~ ^-^ ^ 10^(1 + z). At z = 3 this yields a spectral resolution of 50 times higher 
than the designed value of the spectrograph aboard FUSE. 

The evolution of the mean transmitted flux is typically represented by the effective optical 
depth, defined as Tg// = —ln{F). As discussed in Paschos & Norman (2005) (PN), the effective 
optical depth is biased by high transmission gaps in the pixel fiux distribution and therefore will 
systematically yield lower values when compared to the mean optical depth. The latter is defined 
as the raw average of the pixel optical depth per redshift interval, Tmean = Npx-nios ^T°^'^i ^^'^^i' 
where Tij is the pixel optical depth at redshift Zi in the line of sight index LOS = j. In Figure ([TO]) , 
we show the redshift evolution of the the optical depth of the 304 A line using both representations, 
mean (left panel) and effective (right panel), which suggests a rather smooth evolution leading to 
Hell reionization at z ~ 2.5. Error bars show the la standard deviation of the mean flux and 
optical depth values between different lines of sight. In the effective optical depth case, standard 
deviation at more than 100% the mean fiux value at z > 5 yield negative lower bound values. A 
cutoff in the mean fiux of Fmin = W~^ was therefore imposed in order to be able to compute the 
natural logarithm. 

The error bars also indicate a large degree of variance in the data at redshifts z > 3.5 — 4 
consistent with PN which discussed properties of Hi transmission during hydrogen reionization. 
There, the large degree of variance during reionization was attributed to the presence of high 
transmission gaps associated with underdense regions in the IGM which ionize first. In this picture 
of inhomogeneous He ii reionization, the high transmission segments at high redshifts are associated 
with Hem bubbles that the lines of sight intersect as they are cast through the simulation box. 
However, this is not inconsistent with the conclusions reached in PN because the Hem bubbles 
primarily extend in underdense to mean density cosmic regions. 

Estimates of the observed Heii effective optical depth, shown in Figure (jlOp . are comprised 
of the six known and analyzed sightlines todate: HS1700+6416 (Fechner at al. 2006; Davidsen 
et al. 1996), HE2347-4342 (Zheng et al. 2004b; Kriss et al. 2001), HS1157+3143 (Reimers et 
al. 2005), PKS1935-692 (Anderson et al. 1999), Q0302-003 (Heap et al. 2000) and SDSSJ2346- 
0016 (Zheng et al. 2004a). Although very few in number they suggest a lack of Heii Lya forest 
transmission at redshifts above z ~ 3 and the presence of a trough there consistent with the 
hypothesis of late He ii reionization. Upper bound estimates suggest a steeply rising optical depth 
above z ~ 2.8 while lower bound estimates infer a smoother reionization transition. Lower redshift 
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sightlines towards HS1700+6416 (Fechner et al. 2006) and HE2347-4342 (Zheng et al. 2004b; 
Kriss et al. 2001) suggest an optical depth of about 1 between z ~ 2.3 and z ~ 2.7. In this 
work, we have calculated that between z = 2.6 — 2.4, at z = 2.5, a mean transmitted flux value at 
F = 0.304 lb 0.002 which corresponds to r^-^-^ = 1.190 ib 0.007. The error to the mean values is 
due to sightline-to-sightline variance and is equal to l(T/\/iVios- Observed effective optical depth 
estimates are r^^/ = 0.91 ±0.01 for HE2347-4342 averaged over z ~ 2.3 — 2.7 and a range between 
T^f/ = 0.74 ± 0.34 and r^f/ = 1.06 ± 0.18 at z ~ 2.45 for HS1700+6416 (Fechner et al. 2006). 
If reionization is completed by z ~ 2.5, as suggested by the observed lines of sight, that would 
correspond to an optical depth of T^ff ^ 1. Our calculation is about 1.5<t above that value at 
that redshift. Furthermore, the computed redshift evolution of the optical depth is smoother than 
suggested by observations. We attribute such differences to the limitations of our treatment and 
most notably ignoring diffuse emission due to recombinations to the ground state of He ii . Such 
recombinations increase the ionization of Heii and therefore lower the opacity. The slope of the 
redshift evolution is affected by such omission because of the diffuse emission's dependency on the 
gas temperature, oc T~2 . At earlier redshifts, the temperature is lower and therefore recombinations 
to the ground state may contribute a significant amount of Heii ionizing radiation. Nonetheless, 
due to the lack of a self-consistent calculation we do not know how large that effect would be on 
the slope of the optical depth redshift evolution. We do anticipate that the effect diminishes with 
redshift due to photoheating which raises the IGM temperature. 



From the right panel of Figure (jlOp we can infer a general agreement between the observed and 
computed values at redshifts that evidently sample the tail of the He ii reionization epoch. That 
leads to the conclusion that the final opacity distribution in He ii is largely insensitive to the details 
of our numerical setup. Direct photo-ionizations due to a rising number density of QSO sources 
appears to be adequate to yield a value of the He ii abundance close to the observed one at late 
redshifts. 



6.3. Line Statistics 

The largely ionized Heii by z=2.5 in our computational volume and the update in the gas 
temperature due to photo-heating, allows for the derivation of two standard statistical properties 
of a Lya forest, the column density and b-parameter distribution. Our grid resolution is too coarse 
to draw any absolute conclusions about these distributions. In order to resolve the Hi Lya forest, 
a grid resolution of ~ 40 kpc is required (Bryan et al. 1999). A larger grid resolution may not be 
required for He ii lines if they primarily form in underdense regions. However, comparison with a 
resolved H i forest is desirable. Our intent in this work is to look at the relative differences between 
the two forests and draw conclusions that may stand the test of an improved grid resolution in a 
future simulation. 

Between z=2.6-2.4, at z = 2.5, we identify the Heii and Hi Lya lines and compute their 
column density and b-parameter width. The line identification was performed using the method 
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described in Zhang et al. (1998). We count lines per constant logarithmic column density bins 
{ANx = 0.25) and normalize per column density and redshift interval, dz = 0.2, for the two 
species, ^^^j^ (X = Hi , Heii ). The result is shown in the left panel of Figure pT]) . The 
Hi distribution flattens at LogNfji ~ 12.5 which is due to the coarse grid resolution. A similar 
flattening occurs at LogNueii ~ 13.5. In addition, the Hen column density distribution shows a 
decrease in the slope and large error bars as the column density increases. The error bars refer to 
the la standard deviation of the identified line counts per logarithmic bin. By LogNneii ^ 17.5 
the dramatic increase in the error and hardening of the slope is interpreted as the effect of regions 
of Hell that are not fully ionized by z=2.5. Such regions occupy about 10% of the computational 
volume and are a source of large optical depths and therefore column densities in lines of sight that 
pass throughO them. The two profiles suggest that for column densities in the range 13.5 - 15.5 
we identify 10-100 times more Heii than Hi lines. This result is consistent with a high resolution 
cosmological simulation that resolves the two forests (Zhang et al. 1997). 

Power law fits to the column density distributions in Figure (|11|) . fx{N)dNx = f3N^'^-^dNx, 
yield aui = 1-89 it 0.14 and aneii = 1-41 ± 0.12 in the logarithmic column density range of 13.5 - 
15.5. The error estimates are derived from the least squares fit and do not take into consideration 
any propagated error in the individual column density bins which is sensitive to the bin size. The 
Lya forest is well studied in the literature, however results that are pertinent to the slope of the 
column density distribution typically refer to either the range of logN}{i{cm~'^) = 12.5 — 14.5 
or to an average slope between the minimum to maximum value in the column density sample. 
From distributions obtained in simulations (Zhang et al. 1997; Jena et al. 2005) and observations 
(Petitjean et al. 1993; Ranch et al. 1997; Kirkman et al. 1997) we have estimated a range 
in the slope of the distribution between Nhi = 10^^'^ — 10^^'^ cm~^ of aui = 1-7 — 1.85 for 
redshifts between z=2-3. The slope for the Hi column density distribution obtained in this work is 
slightly above. That is primarily due to the poor grid resolution that underestimates the number 
of Hi absorbers with high column densities. To a lesser extent, we also expect a steepening in 
the distribution because we overestimate the H i opacity for lower column density absorbers in our 
numerical setup. However, we note that within the error estimate of the power law fit, the value 
of the slope we computed here overlaps with the range of published values. 

The estimated slope of the Heii Lya column density distribution, in the optical thin limit, 
using uniform photoionization models, ranges between aneii ~ 1-5 — 1.6 in the redshift interval 
z=2-3 and for logNneii = 13.5 — 15.5 (Zhang et al. 1997). In this work, we have calculated the slope 
of the distribution at z = 2.5 to be aneii = 1-41 which is below the previous estimates. As we shall 
show below the Hell absorbers correspond to physically extended IGM structures and therefore 
the Hell column density distribution is less sensitive to our coarse grid resolution than hydrogen. 
The smaller slope is due to the opacity effects of the inhomogeneous ionizing radiation. Overdense 
regions will have higher opacity to Heii ionizing radiation compared to a uniform photoionization 
model. That in turn results in lines migrating to higher column density bins which softens the 
slope of the distribution compared to an optically thin calculation. 
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We now proceed to investigate the effects of non-uniform photoheating in the distribution of 
the hue broadening widths which is the observable that closely traces the thermal state of the IGM. 
The b-parameter distribution, computed here as the fraction of lines per km/s is shown for the two 
species on the right panel of Figure (jlip (solid lines). The distributions were computed from lines 
with column densities in the range LogNx = 13.5 — 15.5 for both species. The selected range aims to 
avoid the undersampling of weak lines due to the coarse grid resolution. For reference, we also plot 
the b-parameter distribution computed from a uniform photoionization optically thin simulation 
due to a Haardt & Madau (2001) UVB model from a mix of quasar and galaxy populations (dashed 
lines). The median line broadening widths derived in the non-uniform UV case are ^^g^ = 34.23 
kms~^ and h^J^ = 28.16 kms~^ for the Hi and Heii lines respectively. For comparison, the median 
values in the uniform UVB case are i^ed(^) ~ 32.98 kms~^ and h^^J/{u) = 26.89 kms~^. When we 
compare the non-uniform and uniform results we note that the peaks of the distributions are offset 
by ~ 1.25 km/s in both species. This is expected according to the right panel of Figure ^ where 
the median temperature per logarithmic overdensity bin is plotted. Differences in the overdensity 
dependent median temperature distribution arise at overdensities A < 3. If we adopt a relation 
between Hi column density and overdensity in the form of A > 10(yj#5-)^'^(— |^)~^ (Schaye et al. 
2003) then we can compute that at z=2.5 a Hi column density of 10^^'^ cm~^ is due to overdensities 
of A > 1.5. At such overdensities we predict the biggest difference in the temperature between the 
non-uniform versus the uniform cases. The increase in temperature by a factor ~ 1.6 is consistent 
with a shift in the peak of the b-parameter distribution. Note however, that the distributions 
at broadening widths larger than 50 km/s are indistinguishable which may be due to collisional 
ionization dominating the line opacity. 

Since the forests are unresolved on the simulation grid, we make no claim on whether the 
computed b-parameter values have any relevance to the real universe. We can however compare the 
median values of the two species. Relative to hydrogen the median values infer that on average the 
Hell lines have about ~ 82% the broadening width of the Hi lines. For pure thermal broadening, 
the heavier by a factor of 4 helium atoms would yield a width only half the corresponding size of 
the hydrogen line, b^ = 2^^ ■ Our calculation suggests that the broadening due to the Hubble 
differential flow and peculiar velocities dominate the line formation in the Heii forest. That in turn 
would indicate that the Heii lines primarily form in underdense and cosmic mean density regions. 
Our result is consistent with the conclusions in Zheng et al. (2004) where they calculated that 
along the HE2347-4342 transmission line bneii — O.dbbni which also supports a Hubble dominated 
absorption line broadening for He ii . 



6.4. ?7-Parameter Evolution 

We conclude this section, by computing the redshift evolution of the ry-parameter through 
the R- factor, which are defined in ^2.2i In Figure (J12p . we plot the redshift evolution of the rj- 
parameter in the range 3.5 < z < 2.5 in redshift intervals of 6z = 0.2, averaged over all lines of 
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sight. Direct computation of the quantity not only requires knowledge of the column density for 
each individual line but also knowledge of the local association between the H i and He ii absorption 
features. However, because the transmission of the Heii Lya line exhibits a trough at redshifts 
2; > 3, obtaining the column density value for individual He 11 lines is not feasible there due to 
the inability to deblend the absorption features in the spectrum. In addition, whatever correlation 
between the H i and He 11 lines exists at z < 3 it would be incomplete because not all H i lines can 
be associated with He 11 absorption features and vice versa (Kriss et al. 2001). Our biggest problem 
though in this calculation lies in the separate treatment of hydrogen and helium ionization. The 
latter coupled to a low resolution simulation can be a source of significant bias against the true 
correlation between the He 11 and Hi absorbers. 

Therefore, we will approximate here the line of sight and redshift interval average of the rj- 
parameter as < log{ri) >~< /og(4 x ^^"^^ ) >=< log{4: x R) > where R stands for the R-factor. 
Through this approximation, the identification of individual lines and the derivation of column 
density through their gaussian profiles are not required. Instead, in each redshift interval we 
compute the average of the ratio between the local He 11 and H i pixel optical depths along the 
transmission line. The values for each line of sight are then averaged to obtain a mean value for the 
ry-parameter. We plot the redshift evolution of such calculation in Figure (fT2]) . Because the errors 
due to pixel averaging are large, they are ignored. Instead, the error bars in the figure represent the 
lo" standard deviation due to line of sight averaging. For reference, we overplot the computed value 
of the 7y-parameter in the spectrum of HE2347-4342 from Kriss et al. (2001) and also data from 
Fechner et al. (2006) (blue circle) towards HS1700+6416. We obtained the latter by computing 
the mean logr] from the published data in the range z=2.3-2.75. At z=2.5 our computed value 
logrj = 2.1 lb 0.1 is above the observed values which is mainly due to our larger estimate of the 
Hell Lya optical depth at that redshift. As with the optical depth calculation, the limitations of 
our treatment result in an overestimate. However, we came close enough to the observed values 
at redshifts later than the Hen transmission trough, to reinforce our conclusion that our limited 
calculation is able to reproduce the much of the observed signatures of He 11 reionization. 



7. Summary & Conclusions 

We have simulated the late inhomogeneous reionization of He 11 by quasars and the attendant 
photoheating of the IGM including opacity effects. We post-process baryonic density fields from 
a standard optically thin IGM simulation with a homogeneous galaxy-dominated UV background 
which reionizes Hi and Hei at z=6.5. Quasars are introduced as point sources throughout the 
100 Mpc simulation volume located at density peaks consistent with the Pei luminosity function. 
We assume an intrinsic quasar spectrum J{i^) oc i/^^-^ and a luminosity proportional to the halo 
mass. We evolve the spatial distribution of the He 11 ionizing radiation field at hz/ = 4, 8, and 
16 Ryd using a time-implicit variable tensor Eddington factor radiative transfer scheme which we 
describe. Simultaneously, we also solve for the local ionization of He 11 to He iii and the associated 
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photoheating of the gas. This distinguishes our calculation from the work of Sokasian, Abel and 
Hernquist (2002) who simulated inhomogeneous Hell reionization but did not study the thermal 
evolution of the gas. 

We find that the cosmic evolution of the QSO population causes the individual Hem regions to 
overlap and subsequently ionize 90 % of the volume at He iii ionization fractions tp^ = ^^"' > 10~^ 
by z ~ 2.5. In addition, to the Heii and Hem number density calculation, we also update the local 
gas temperature due to the thermalization of the photoelectrons ejected by the HeH ionization 
process. As expected, the temperature is higher compared to the unprocessed simulation with 
the difference increasing rapidly at redshifts ^; ^ 3. Relative to a self-consistent optically thin 
simulation where He ll is also photoionized by a homogeneous UV background, we find that optical 
depth effects result in an increase in the temperature of the intergalactic medium at the cosmic 
mean density level by a factor of ~ 1.7 at z = 2.5. The results of our temperature calculation are 
consistent with analytic and numerical predictions of the HeH heating effect (Abel & Haehnelt 1999; 
Schaye et al. 2000; McDonald et al. 2000). Finally, we trace the redshift evolution of the Heii Lya 
transmission using randomly casted synthetic spectra through the simulated volume. The analysis 
of the mean transmission allows for the derivation of the effective optical depth of the 304 A line. 

We have calculated that at z = 2.5ib0.1, the average of pixel flux among all lines of sight yields 
a value for the mean transmission oi F = 0.304 it 0.002 which corresponds to t^I = 1.190 ±0.007. 
The length of each sightline (Az = 0.2) is longer than the size of the simulation in order to 
minimize the transmission variance across the redshift path. We compute the error only due to the 
sightline to sightline variance and find that at z=2.5 the la standard deviation is 11% the mean 
Lya flux. When we compare to estimates from Heii forest spectra observed with the FUSE (Far 
Ultraviolet Space Explorer) satellite and HST at the same redshift interval we find that the our 
value of the effective depth is comparable but slightly above the observed values. We attribute the 
disagreement to our approximate treatment of the inhomogeneous ionizing radiation that ignores 
the diffuse component and only focuses on the point sources' input. 

In addition to the mean transmission estimate, we also compute the column density and b- 
parameter distribution of the identified Hen Lya lines and compare them to the corresponding 
statistical properties of the HI Lya lines. We find that an optically thick calculation results in 
extra heating that shifts the b-parameter distributions to higher broadening widths. In addition 
and within the limitations of our coarse numerical resolution, which does not resolve the Hi forest, 
our calculation shows that the median Heii b-parameter at z = 2.5 is 82% the HI broadening 
width. The latter suggests that the Hubble differential flow may be the dominant line broadening 
mechanism which in turn indicates that Hell lines primarily form in underdense to cosmic mean 
density regions. 
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A. Species Concentration Equation in an Expanding Universe 

For a single species medium and in proper coordinates the rate equation that governs the 
abundance of ionization state i is given by the following equation: 

hi = -S-rii + T,jPijnj (Al) 

In Equation (jAip . rij denotes the number density of ionization states of the species in ionization 
or recombination coupling to state i. The ratio - is a function of redshift and is equal to the redshift 
value of the Hubble constant. 

Dividing by the total number density of the species, nx, we get: 



hi/riT = —3—ni/nT + T^jPijUj/riT ^ 

hi/riT = -3-yi + T,j(3ijyj (A2) 

The LHS of Equation (|A2p . can be rewritten as follows: 



i = yi- S-Ui (A3) 

a 

Plugging the result back into Equation (|A2|) . we get in = T^jPijTjj, which shows that when 
solving the rate equation using ionization fractions, instead of number densities, the equation 
is independent of the Hubble expansion term. However, coefficients (5ij that are pertinent to 
recombinations to state i depend on the electron density which is computed on the proper frame 
of reference. 



B. Effect of Late Reheating on Hydrogen Ionization Balance 

In this calculation, we estimate the neutral hydrogen fraction if an additional source of ion- 
ization and heating is introduced due to the full photoionization of He ii by a local UV radiation 



1 


d 




^ f 


rii , 




riT 


m""' 


d 


¥t^- 


1 


— n. 






dt 


yi + 





33 



followed and the ejection of an additional electron per helium atom, assumed to be already singly 
ionized. The unmodified neutral fraction is labeled as Case (I). Case (II) represents the updated 
hydrogen abundance. 

In Figure ()13p . Equation ([6]) was solved for two values of temperature increase, T jT = 
1.5 — 2. On the left panel we show the neutral hydrogen fractions in Case (II) for a range of 
Case (I) inputs. The open squares (circles) represent the smallest (largest) temperature increase. 
The solid line stands for no difference between the two cases (slope of 1). The calculation shows 
that the largest temperature increase results in a greater departure from the straight line, although 
the degree of such departure becomes smaller at larger neutral fractions. 

This is shown on the right panel of Figure (I13p where we plot the percentage change in the 
hydrogen neutral fraction. We see a fixed degree of change between x{ = 10^^ — 10^^'^ and a rapid 
decrease at xi ^ 10^^'^. For the range of temperature increase used in the calculation we can then 
determine that a total He ii photoionization would lead to an adjustment of the hydrogen neutral 
fraction by 12-24%. This corresponds to an overestimate of xi, in the absence of Heii photoheating, 
by « 14 - 32%. 

The degree at which any temperature change affects the hydrogen ionization balance is sensitive 
to the functional form of the recombination coefficient. In general, the coefficient steepens at 
temperatures above T ~ 10^ K where the slope becomes a function of temperature itself. The 
results shown in Figure (J13p are derived with a constant slope of /? = 0.5 which is roughly valid at 
T ~ 10^ K. To explore the effects of the slope on the adjustment of the hydrogen neutral fraction, 
due to He ii photoheating, we show on the left panel of Figure (J14p , the results of the calculation 
for an input neutral fraction of Xi = 10~^, /? = (0.4,0.5,0.6) and for a range of temperature ratios 
(0.5-2.0). 

The upper dashed line at T /T > 1 represents the largest /? value (0.6) and indicates that 
a steepening of the temperature exponent leads to an increase in the neutral hydrogen fraction 
adjustment. This trend is reversed at T^^ /T^ < 1 (cooling). The latter temperature adjustment 
range is not related to any Heii photoheating eff'ects but it was a numerical investigation of the 
solution. However, one can make the argument, depending on whether the UVB drops after a 
certain redshift, that additional cooling terms, such as collisional ionization cooling {Hejj + e~ ^ 
Hem + 2e~) and recombination radiation cooling {Hem + e~ — > Heji + 7), might lower the 
temperature below the unperturbed value. However, such possibility is small at least up to z ~ 2. 

We also note from the left panel of Figure (J14p that not all values of a temperature increase 
lead to a positive neutral hydrogen fraction adjustment (1 V > 0). In the f3 = 0.5 case 

Xi 

(solid curve) a small increase of the temperature by less than 15adjusted to a larger value and 

// 

therefore 1 V < 0. This a result of the increased electron number density due to the additional 

xi - 
photoelectron ejected from the He 11 atom. In this narrow range of temperature increase the greater 

electron number density dominates the recombination rate, which consequently increases to shift 

the ionization balance towards more neutral hydrogen. 
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This effect depends on the initial neutral fraction, as shown on the right panel of Figure ([T] 
There, we plot the family of solutions of Equation dS]) (/3 = 0.5) for the range of neutral fractions 
(10~^ < xi ^ 10~^) and temperature adjustment factors (0.5 < T^^ /T^ < 2.0). The solid line 
represents the solution obtained for x{ = 10~^ from the left panel. The plot shows that as the 
neutral fraction increases the adjustment decreases for T^^ /T^ > 1. As discussed in the previous 
paragraph, the range of T /T where the electron density dominates the recombination rate shifts 
to the right towards lower ionization fractions. 



C. Radiative and CoUisional Rates 

In this section we document the radiative and collisional rates we used in our simulation. For 
the full range of chemical reaction rates incorporated into the cosmological hydro code Enzo please 
see Anninos et al. (1997) and Abel et al. (1997). 



C.l. Collisional Ionization and Radiative Recombination of Singly Ionized Helium 

Fits are accurate within 1% in the temperature range 1 — 10^ K. 

He+ + e- ^ He++ + 26" Abel et al. (1997) 

In the following fit temperature T is in units of eV: T = hqqq p^ 

TY'/ne = h = exp[-68.71050990 + 43.93347633/n(r) 

-18.4806699/n(r)2 + 4.70162649Zn(T)^ - 0.76924663/n(r)^ + 

8.11042 X 10"^/n(T)^ - 5.32402063 x 10"^/n(r)'^ + 

1.97570531 X 10"'^;n(r)'^ - 3.16558106 x 10"^/n(r)^] cm^s~^ 

He++ + e- ^ He+ + 7 Cen (1992) 

a%^h = 3.36 X 10-i°r-V2(_r_)-o.2(i + (^)0.7)-i amh-^ 

Temperature is expressed in K. 

C.2. Photoionization cross section of Hydrogen and Helium II 

He+ + 7 ^ He++ + e~ 
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The cross section for both rates are expressed through the same functional form because singly 
ionized helium is a hydrogen like atom (from Osterbrock 1974). 

A ( u \4.exp[4—A(arctane)/t] 

'^ ~ 'Z^^'u^' l-exp{-2-n/e) 

A = 6.30 X 10-1^ cm^, e = (^ - 1)^/^ and uz = 13.6 x Z^ (Z = 1, 2). 

For consistency with the notation used in the main text, we need to clarify that we we label 

For numerical convenience we use an approximation to the full Hell cross section in the form 
(^Heii = j(= '^°Heii)ihu^)~^i where /if3 = 54.4 eV. Such an expression allows an analytic calculation 
and fast reconstruction, via a polynomial fit, of the energy density profile between four frequency 
points. The power law approximation of the cross section is estimated to be accurate to within 
K, 3% (Haehnelt k Abel 1999). 
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Fig. 1. — Left Panel: Volumetric rendering of the logarithm of baryon overdensity at z=2.6. Right 
Panel: Evolution in the number of QSO sources versus redshift in this simulation 
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Fig. 2. — Volumetric renderings of the Hem distribution mapped by the Lognneiii at z=3 and 
z=2.6. The ionization cutoff in this visuahzation is set to tp^ = "■'^•^'^ = 10^^ (dark regions). 
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Fig. 3. — Left Panel: Slice of the Hem (darkened regions) and Hen 3D mass distribution at z=2.6. 
At that redshift, the volume filling fraction of the Hem is VFF ~ 0.90. Right Panel: Redshift 
Evolution of the VFF for XHeiii ^ 10^"'. The error bars result from the uncertainty in the location 
of the cumulative I-front due to finite grid cell size. 
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Fig. 4. — Left Panel: Redshift evolution of the volume averaged mass fraction of singly helium. 
Overplotted is the corresponding neutral hydrogen mass fraction. Right Panel: The scatter plot 
between Heii fraction and local gas overdensity at z=2.5 yields the mean and median value of 
Log f Hell per logarithmic bin of gas overdensity. Solid (dashed) curves refer to the inhomogeneous 
(homogeneous) case. The median profiles correspond to the curves that peak at Log6 ~ 1 in each 
case (the mean profiles peak at a lower overdensity). 
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Fig. 5. — Left Panel: Redshift of the temperature evolution. Shown as a solid (dashed) line is 
the evolution in the processed (pure) case of the mean temperature in grid cells with overdensities 
between 6 = 1 — 10. Right Panel: The median temperature distribution in each logarithmic 
overdensity interval at z=2.5. It shows that the processed (solid) and unprocessed (dashed) cases. 
The curves show that He ii reionization primarily increases the temperature of the IGM in low to 
mildly over dense range. 
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Fig. 6. — Left Panel: From temperature-density plot in Figure ([5]) we compute the ratio of the 
postprocessed result over the input simulation (dashed) and the optically thin simulation (dashed- 
dot). Right Panel: We plot the distribution of the logarithmic slope -^j^ = 7 — 1 versus the gas 
overdensity density. 
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Fig. 7. — Hell and Hi Lya spectra between z = 2.6 — 2.4, at z = 2.5, along a randomly casted 
sightline through the computational volume. The upper panels refer to the HI Lya transmission. 
The lower panels refer to the Hell Lya transmission. The horizontal axis is converted to observed 
wavelength A = Xoil + z), where Aq is the restframe wavelength of the resonant Lya scatter. 
Left Panels: Uniform UVB ionizes and heats the IGM in the optically thin approximation. HeH 
ionization and heating is included in this calculation. Right Panels: The uniform UVB only ionizes 
Hi to Hii and Hei to Heii . Photoheating therefore due to the uniform UVB only refers to the 
above processes. The non-uniform point source radiative input from QSO type sources is used to 
photoionize He ii to He iii and to calculate the subsequent photoheating. Although the H i spectra 
seem similar between the two sets of calculations, ignoring the feedback on Hi abundance by the 
Hell ionization results in underestimating the mean Hi Lya transmission, compared to the full 
simulation, by « 8%. 
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Fig. 8. — Impact parameters of the random sightlines to the point sources in the box at z = 2.5. At 
each point along the path of the sighthne, we compute the distance to the closest point source. As 
the trajectory crosses the computational volume, the closest source changes depending on position. 
On the lower left part of the figure, we plot this proximity distance for a sightline that comes the 
closest to a point source, about one grid resolution element. On the top left panel, we show the 
He II transmitted flux for that sightline. Crosses mark the locations along the path that are closest 
to the point source. The alternating proximity sources are indicated by the alternating parabolic 
profiles. On the right panel, we show the probability distribution of the minimum distance to the 
sources, which we define as the impact parameter, from all random sightlines. The asymmetry in 
the distribution is small, which is due to the near isotropic distribution of high density peaks at 
the scale of the simulated volume. 
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Fig. 9. — On the left panel, we show the scatter plot between the average impact parameter for 
each sightline in respect to the distribution of the ionizing point sources versus the sightline-specific 
mean flux. The red line is power law fit to the data with a slope of s = —0.16 it 0.04. The blue solid 
(dashed) line shows the average {2a) Heii Lya transmission from all sightlines at z = 2.5. The 
small negative correlation r = —0.24 suggests that the mean flux of a sightline is largely insensitive 
to the average proximity in respect to the distribution of ionizing sources along it's path. On the 
right panel, we show the locally absorbed flux (y-axis) at the point of closest proximity to the 
ionizing source shown on the x-axis. Individual points correspond to the absorbed flux at the pixel 
of closest distance. The red (blue) histogram is the median (mean) absorbed flux per constant 
logarithmic bins. Solid and dashed blue horizontal lines are the level of the mean absorption in 
our sightline sample and the 2a standard deviation respectively. The solid line is a Fermi-Dirac 
function that fits the median (red) histogram for impact parameters less that 10 Mpc comoving. 
The half point of the Fermi-Dirac function is at do = 3.1 Mpc comoving (shown as a vertical black 
solid line) and it has a skin width of 0.63 Mpc comoving. Impact parameters along the sightline 
of < 4.5 Mpc comoving, the point where the absorbed flux falls below the 2a level of the mean 
absorption, show a steep sensitivity in the absorbed flux due to the quasar proximity. 
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Fig. 10. — We compute the mean optical depth and effective optical depth at selected mean 
redshift points, from z = 6 — 2.5 every 0.25 redshift units and in intervals of 6z = 0.25 about the 
mean redshift. Left Panel: Redshift evolution of the mean pixel Lya Heii optical depth averaged 
over all lines of sight. The error bars are la standard deviation to the sightline average optical 
depth. Right Panel: Redshift evolution of the Lya Heii effective optical depth. The error bars 
in the simulation results are sightline-to-sightline errors estimated from the standard deviation to 
the average transmission from sightline-specific mean flux. Overplotted are sightline measurements 
along HE2347-4342 (red triangles), HS1700+6416 (blue circles), HS1157+3143 (blue diamond), 
PKS 1935-692 (blue cross), Q0302-003 (green square) and SDSSJ2346-0016 (orange star). Error 
bars to these points refer to error estimates of the pixel flux along the sightline specific to the 
quasar. The HE2347-4342 data from Zheng et al. (2004b) have been resampled to a larger bin 
size to make it comparable to our data. The redshift profiles of the optical depth point statistics 
suggest a smooth evolution of the He ii reionization epoch. 
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Fig. 11. — Left Panel: Column Density Distributions, f{N) — ^^^j^ 
between z = 2.6 and z = 2.4. Circle and square points correspond to the Hi and Hen column 
density distribution respectively. Solid (dashed) straight lines (in log-log) show the power law fit 
to the Hi (Heii ) distribution between column densities Nx = 10^^'^ — 10^^'^ cm~'^. The error 
bars show the Icr standard deviation to the number of lines per logarithmic column density bin 
equal to 0.25. Right Panel: Broadening width distributions at z = 2.5 for Hen and Hi . The 
distributions are shown as the fraction of absorption lines per km/s in bins of 2 km/s. Black and 
red colors correspond to H i and He n distributions respectively. The dashed lines show a reference 
self-consistent calculation in the optical thin limit using an ionizing uniform UVB. The peaks of 
non-uniform calculations are offset by ~ 1.25 km/s to the right of the uniform results. 
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Fig. 12. — Redshift evolution of the r] parameter in this simulation. Results are shown as < logrj > 
per redshift interval and are obtained as the line column density ratio of Heii over Hi at low z 
{z < 3) and estimated from the ratio of optical depths at higher redshifts. The filled triangle shows 
the average estimate between z = 2.3 — 2.7 towards HE2343-4342 from Kriss et al. (2001). We 
also show the data from Fechner et al. (2006) (blue towards HS1700+6416 collapsed into a single 
redshift interval from z=2.3-2.75. Our calculation yields systematically larger values, a result due 
to the larger estimate of the effective optical depth hy z = 2.5. 
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Fig. 13. — Left Panel: Logx{ (input) vs. Logxi (output) for T /T = 1.5 (squares) and T /T 
2.0 (circles). Right Panel: Percentage change in the hydrogen neutral fraction vs. Logx ■ 
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Fig. 14. — Left Panel: Percentage change vs. T /T for three values of the exponent in a^^ oc T^'^. 
Solid line: [5 = 0.5. Upper dashed line at T /T > 1: /3 = 0.6. Lower dashed line at T jT > 1: 
P = 0.4. Right Panel: Percentage change vs. T /T for 10^^ < Xi ^ 10~^. Curves on the upper 
quadrant from upper to lower (dashed lines) represent an increase in the hydrogen neutral fraction. 



